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Users  Guide  for 
ESD  LORAN  Grid  Prediction  Program 


1.  INTRODl CTION 

The  ability  to  precisely  deliver  ordnance  or  men  and  material  under  all 
weather  conditions  in  an  adverse  battlefield  environment,  is  one  of  the  most 
severe  requirements  imposed  on  U.  S.  Worldwide  Tactical  Air  Forces.  A key 
element  in  meeting  this  requirement  involves  obtaining  and  providing  position 
data  in  a specified  coordinate  system,  to  allow  navigation  of  tactical  aircraft  to 
desired  locations  with  sufficient  accuracy  for  rendezvous  or  weapon  release.  The 
LORAN  radio  navigation  system  is  being  relied  upon  more  and  more  as  a principle 
source  of  navigation  information  in  tactical  airborne  operations.  Position  indica- 
tion is  given  in  terms  of  LORAN  TD's  (time  difference)  which,  because  of  propa- 
gation anomalies,  do  not  correspond  precisely  to  earth  fixed  geodetic  coordinates. 
Therefore,  each  LORAN  chain  requires  a grid  prediction  for  its  coverage  area. 
Such  a grid  prediction  computer  program  package  has  been  developed  at  RADC  ' 
EEP  and  is  described  herewith.  This  manual  contains  sufficient  information  to 
enable  the  experienced  programmer  to  understand  the  programming  aspects  of 
LORAN  grid  prediction  and  includes  a detailed  functional  description  and  its 
operation. 


(Received  for  publication  8 December  1977) 
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The  operation  of  the  entire  l^ORAN  Grid  Prediction  System  from  obtaining  the 
required  input  data  off  maps  to  the  calculation  of  the  output  time  difference  is 
illustrated  in  Figure  1.  The  system  has  been  divided  into  the  four  following  parts: 

PART  I:  Formation  of  map  digitized. 

PART  II:  Calculation  of  surface  impedance  data  base. 

PART  III:  Field  LORAN  prediction  package. 

PART  IV:  Updating  data  base  by  comparison  of  calculated  and  experimental 
results. 

A sufficient  description  of  the  various  computer  program  is  presented  so  that 
the  user  can  obtain  the  required  LORAN  coordinates.  The  data  base  is  generated 
in  PARTS  I and  II  and  updated  when  measured  data  is  available  in  PART  IV.  In 
PART  III  of  the  systems,  the  main  set  of  calculations  are  performed  and  this  is 
described  first  in  Section  3.  The  technique  for  translating  the  earth's  electrical 
properties  into  a surface  impedance  and  properly  sequenced  onto  a rapid  access 
magnetic  disc  is  described  in  Section  4.  Brief  description  of  map  data  digitization 
and  system  tuning  or  data  base  updating  is  given  in  Sections  5 and  6,  respectively. 


PART  1 


PART  11 


PART  111 


SOIL 

GEOLOGY 

ELEVATION 


DATA 


Figure  1.  Total  LORAN  Grid  Prediction  System 
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.1.  DKSr.KlITION  OK  I*  MIT  III 
3.1  (>ver\ie» 

A block  diagram  of  the  PART  III  LORAN  prediction  package  id  shown  in 
Figure  2.  Its  purpose  is  to  furnish  the  LORAN  TD's  for  a desired  target  given  the 
ground  properties  of  the  system  coverage  area.  For  each  LORAN  chain,  the  geo- 
detic location  of  the  master  and  two  slave  transmitters  must  be  known  in  addition 
to  the  two  slave  emission  delays.  A magnetic  tape  of  the  ground  electrical  prop- 
erties for  the  given  coverage  area  is  also  required.  The  latitude  and  longitude  of 
the  desired  target  and  delivery  altitude  are  inserted  into  the  program  and  a time 
of  arrival  (TOA)  from  each  of  the  transmitters  is  computed.  Subtracting  the 
master  TOA  from  each  of  the  slave  TOA's  yields  two  TD's  which  determine  the 
LORAN  coordinates  of  the  target. 


TAPr  INPUT 
OP  CROUND 
PROPEKTiES 


Figure  2.  LORAN  Prediction  Package  — Part  III 


The  operation  of  this  program  can  be  followed  by  referring  to  Figure  2.  The 
input  tape  contains  the  ground  electrical  properties  which  consist  of  elevation  and 
complex  surface  impedance  for  the  area  of  interest.  This  is  in  the  form  of  a 
matrix  of  data  points  every  30  arc  sec  in  latitude  and  longitude.  The  information 
is  recorded  onto  a disc  for  rapid  access  of  data  between  any  two  points  in  the  serv- 
ice area  as  required  by  the  geodesic  retrieval. 
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The  input  data  cards  contain  the  following:  (a)  transmitter  coordinates  with 
associated  emission  delays,  and  (b)  target  coordinates,  delivery  altitude,  inte- 
gration step  size,  and  type  of  receiving  antenna  for  each  target.  This  data  is  first 
used  in  the  geodetic  calculation  to  determine  path  length  (D)  between  each  of  the 
transmitters  and  the  target,  and  to  determine  points  along  the  geodesic  path 
governed  by  the  distance  increment  or  integration  step  size.  The  points  along  the 
geodetic  path  are  input  to  the  geodesic  retrieval,  which  in  turn  ot  tains  the  ground 
electrical  properties  of  the  points  from  the  matrix  on  the  disc.  This  information 
is  used  to  calculate  the  time  correction  or  secondary  phase  factor  due  to  the 
‘ decrease  in  propagation  velocily,  compared  to  free  space,  when  a signal  propagates 

. over  the  earth's  surface.  Time  of  arrival  calculations  can  then  be  made  from  the 

following  equation: 

TOA  = - D I t + ED  (1) 

c c 

j where 

n = atmospheric  index  of  refraction  - 1.000338. 

c - velocity  of  light  2.  997925  X 10^  meters/second. 

I D = length  of  geodetic  path  from  transmitter  to  receiver  in  meters. 

j t = time  correction  or  secondary  phase  factor  for  propagation  over  a given 

^ path  length  in  psec. 

I ED  = emission  delay  in  psec. 

i 

: Three  such  calculations,  one  from  each  transmitter  to  target,  ai-e  required 

for  each  prediction.  Subtracting  the  master  TOA  from  each  of  the  slave  TOA's, 
one  obtains  the  LORAN  coordinate  prediction. 

3.2  Pro^irain  How  C-liarl 

A flow  chart  illustrating  the  operation  of  the  LORAN  Grid  Prediction  package 
is  illustrated  in  Figure  3 and  subroutine  relationships  to  the  driving  program 
LORANCO  is  shown  in  Figure  4.  This  package  consists  of  a deck  with  approxi- 
mately 1300  cards  and  requires  a core  memory  of  120K  base  eight  (8).  The  three 
transmitter  (M,  SI,  S2)  coordinates  with  corresponding  emission  dela\'S  are  read 
into  the  program  from  the  input  data  deck.  For  computation  purposes,  the  geo- 
I graphic  coordinates  are  converted  into  radians  by  a call  to  subroutine  COR  RAD. 

The  target  coordinates  are  then  read  into  the  program  w’ith  corresponding  informa- 
^ tion  on  delivery  altitude  in  meters,  computation  step  size  in  kilometers,  and  t>'pe 

of  receiving  antenna  (Electric  dipole  = 1,  Loop  = 0).  Similarly,  these  coordinates 
are  converted  to  radians.  The  program  is  terminated  when  the  step  size  (ADELS) 
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READ  THREE  Xmtr  Coordinates  and  Coding  Delays 


READ  Receiver  Coordinates,  Altitude,  Step  Size  and  Antenna  Type 

(Z  Alt)  (ADELS)  (NUT) 


ADELS, LE.O  GO  TO  700  EXIT 


CONVERT  LAT/LONG  TO  RADIANS 
SUB  CORRAD 


Define  Total  Distance  from  Transmitter  to  Receiver  D 
by  Sadano  SUB  GEODI 


Define  incremental  distance  by  Hufford 
SUB  GEOPTS 


Obtain  Geodetic  Data  from  Disc  Data  Base 
SUB  GEORET 


Calculate  Secondary  Phase  Factor  t^ 
SUB  INEQ 


For  Ground  or  Altitude  Define 


Define  True  Time  Delays  and  PRINT 


Go  To  5 

New  Receiver  Coordinates 


Figure  3.  LORAN  Program  Flow  Chart 


rr- 

I 


ATI 

AT2 


Figure  4.  ESD  LORAN  Grid  Prediction  Program 


is  equal  to  or  less  than  zero.  Therefore,  a blank  card  after  the  last  desired  tar- 
get card  causes  the  program  to  exit.  With  information  on  transmitter  and  target 
location,  the  program  can  now  determine  the  propagation  path.  The  total  geodesic 

distance  for  primary  wave  delay  is  calculated  in  subroutine  GEODI  which  uses  the 
1 2 

Sodano  formulation.  ’ The  coordinates  of  incremental  points  along  the  propaga- 
tion path  are  defined  in  subroutine  GEOPTS.  The  ground  electrical  properties  for 
these  points  are  obtained  through  subroutine  GEORET  which  in  turn  calls  subrou- 
tine GETFLV.  In  the  latter  subroutine,  the  required  points  are  indexed  to  address 
data  for  the  random  access  disc.  Upon  return  from  the  disc,  the  data  is  unpacked 
and  returned  to  GEORET.  With  a detailed  description  of  the  points  along  the  prop- 
agation paths  available,  one  can  now  calculate  the  secondary  phase  factor  or  time 

correction,  t , due  to  the  disturbing  influence  of  the  earth  using  the  Huffordlnte- 
^ 3 4 5 

gral  formulation  equation.  ’ ’ Ninety-five  percent  of  the  compute  time  is  re- 
quired for  this  calculation  in  subroutine  INEQ.  The  time  of  arrival  is  the  sum  of 
the  primary  wave  travel  time,  secondary  phase  factor,  and  emission  delay.  All 
the  required  information  for  this  calculation  over  a given  path  is  now  available  and 
the  resultant  calculation  is  stored. 


I 

J 

1 


I 

j 

I 


Due  to  the  large  number  of  references  on  this  page,  the  references  will  not  be 
footnoted.  .See  references,  page  .83. 
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The  above  procedure  is  repeated  three  times,  once  from  each  transmitter  to 
receiver  by  looping  back  to  statement  15  after  each  completed  TOA.  With  the 
completion  of  the  TOA  calculations,  the  predicted  time  difference  for  a given  geo- 
graphical coordinate  is  calculated  (T,  T2)  in  the  driving  program  LORANCO  and 
output  printed.  Control  is  then  transferred  back  to  statement  5 where  new  target 
coordinates  are  read  into  the  system,  and  the  entire  process  is  repeated  until 
ADELS  is  made  equal  to  or  less  than  zero. 

Figure  4 is  an  additional  flow  chart  illustrating  the  subroutine  relationships 
to  the  driving  program  LORANCO. 

3.3  Data  Input  and  Output  Setup 

3.3.1  INPUT 

Program  LORANCO  requires  a data  deck  and  data  tape  for  operation  loaded 
as  shown  in  Figure  5.  The  data  deck  supplies  the  geographic  location  of  trans- 
mitters and  targets  and  additional  required  constants  such  as  emission  delay, 
altitude,  step  size,  and  type  of  receiver  antenna.  The  data  tape  supplies  the 
ground  electrical  properties  for  the  entire  service  area  covered  by  the  transmitters. 

The  first  three  cards  in  the  data  deck  describe  the  transmitter  input  data  as 
follows : 


Cols 

Data 

Format 

1-2 

Blank 

3-6 

Alpheric  numeric 

3A8 

identifier  for 

transmitter 

26-42 

Latitude  data 

15,  13,  F7.3,  A1 

26-30 

Latitude,  degrees 

15 

31-33 

Latitude,  minutes 

13 

34-41 

Latitude,  seconds 

F7.  3 

42 

Latitude,  N or  S 

A1 

43-58 

Longitude,  data 

15,  13,  F7.3,  A1 

43-4  7 

Longitude,  degree 

15 

48-50 

Longitude,  minutes 

13 

51-57 

Longitude,  second 

F7.  3 

58 

Longitude,  E or  W' 

A1 

59-78 

Emission  Delay 

F20.  3 

All  degrees,  minutes,  and  seconds  are  right  justified  in  their  respective  field. 
Emission  delay  is  in  units  of  psec.  The  order  of  transmitter  cards  are  master, 
slave  1,  and  slave  2.  The  master  emission  delay  will  always  be  zero.  Slave  1 
emission  delay  is  always  less  than  that  of  slave  2. 
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Figure  5.  Deck  Setup  for  Program  LORANCO 


The  next  set  of  cards 

contain  the  target  input  data.  It  is  read  in  as  one  target 

per  card  so  there  will  be 

as  many  receiver  cards 

as  targets  with  the  following 

format: 

Cols 

Data 

Format 

1 -2 

Blank 

3-26 

Alpheric  numeric 
identification  for 
receiver  identifi- 
cation 

3A8 

26-42 

Latitude  data 

15,  13,  F7.3,  A1 

26-30 

Latitude,  degrees 

Id 

31-33 

Latitude,  minutes 

13 

34-41 

Latitude,  seconds 

F7.3 

42 

Latitude,  N or  S 

A1 

43-58 

Longitude  data 

15,  13,  F7.3,  A1 

43-47 

Longitude,  degrees 

15 

48-50 

Longitude,  minutes 

13 

14 


r 


Cols 

Data 

Format 

51  -57 

Longitude,  seconds 

F7.  3 

58 

Longitude,  E or  VV 

A1 

59-70 

Altitude,  meters 

F12.  1 

71-75 

Step  Size,  kM 

F5.  2 

76-78 

Type  of  Antenna 

13 

0 or  1 


The  aircraft  altitude  at  the  release  point  is  specified  in  meters  and  the  distance 
increment  or  step  size  in  kilometers  typically  0.  5 kilometers.  NUT  defines  the 
type  of  receiver  antenna.  For  a vertical  antenna,  a 1 is  placed  in  column  78,  for 
a loop,  a 0 in  column  78.  A blank  card  or  zero  in  columns  71-75  will  terminate 
the  program. 

The  data  tape  contains  information  on  the  electrical  properties  of  the  ground 
and  covers  the  entire  service  area.  Points  selected  outside  this  area  will  cause 


the  program  to  print  "OUT  OF  ACCEPTABLE  RANGE.  FURTHER  CALCULATIONS 
F'OR  THIS  PATH  HAVE  BEEN  DELETED."  The  supplied  magnetic  tape  contains 
5760  records  of  60  words  each  with  each  word  representing  120  data  points  or  one 
degree  of  latitude.  Each  data  point  is  defined  by  a complex  surface  impedance  and 
an  elevation.  The  current  area  covered  is  66  to  14  degrees  in  longitude  and  48  to 
54  degrees  in  latitude.  This  tape  is  read  into  the  random  access  data  base  disc 
unchanged  in  format  with  the  sequential  record  number  on  the  tape  becoming  the 
random  access  index  array  address  number.  No  operation  on  the  tape  is  required. 


3.3.2  OUTPUT 

Program  LORANCO  produces  the  following  printed  output: 

(a)  For  each  transmitter,  an  echo  printout  of  columns  3-58,  of  the  input  card. 

(b)  For  each  receiver,  an  echo  printout  of  columns  3-58  of  the  input  card. 

(c)  For  each  transmitter  to  receiver  path,  geodesic  path  information,  and  a 
printout  of  the  parameters  used  by  subroutine  INEQ. 

(d)  Results  of  the  ground  wave  time  delay  calculations  in  the  form  of  the  list 
NAME  The  list  NAMI  contains  information  on  time  of  arrivals  (TOA),  time  delays 
(TDl,  TD2),  emission  delays  (EDI,  ED2),  geodesic  distances  (DISTSOD),  primary 
wave  times  (TPW),  and  secondary  phase  factor  times  (TIMDUM). 


3.1  I’rniirain  .''uliroulincs  anil  I' uiR  tioiis 
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LORANCO 

COR RAD 
CORDM3 
GEOPTS 


Reads  inputs,  defines  three  paths,  calls  data  base  and  INEQ, 
calculates  TOA  and  TD. 

Degrees  (l.at.  Long)  to  radians. 

Radians  to  degrees  (Lat,  Long). 

Defines  points  along  geodetic  from  transmitter  to  receiver. 


1 


I 

i 

j 

! 
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IMT 

GETELV 

UNPACK 

GEORET 

INEQ 

INDF 

CNEUKEN 

GEODI 

SETUP 

GROUND 

GANG 

INDEX 

WERE 

OMCOS 

FLEAF 


Calculate  first  and  second  derivative  of  elevation. 

Read  geophysical  data  off  disc. 

Unpack  data  from  disc. 

Returns  geophysical  data  to  driving  program. 

Solution  of  integral  equation  for  secondary  phase  factor. 

Induction  field  calculation  (E  or  H field). 

Interpolation  routine  for  integration. 

Calculation  of  total  geodetic  distance  by  Sodano  and  back  azimuth. 
Constants  for  spheroid  to  be  used. 

Introduces  variable  ground  impedance  and  variable  ground  terrain 
into  integral  equation  formulation  of  the  ground  wave. 

Calculates  argument  of  a complex  number. 

Calculates  index  values  of  data  base  variables  from  LAT /LONG. 
Calculates  error  function. 

Calculates  1 - cos  (X). 

Calculates  ground  wave  attenuation  function  over  flat  ground 
using  flat  earth  theory. 


3.5  Subroutine  Description 

3.5.1  LORANCO 

This  is  the  driving  module  for  the  entire  program.  It  reads  the  input  coordi- 
nates of  both  transmitters  and  receivers,  defines  the  geodesic  path,  receives  data 
base  parameters  and  calculates  the  secondary  phase  factors,  time  of  arrivals,  and 
time  difference.  The  first  three  read  cards,  one  for  each  station,  furnish  the 
transmitter  locations  and  associated  emission  delays.  A call  to  CORRAD  changes 
the  units  of  the  input  data  from  degrees  to  radians.  The  first  target  or  receiver 
geographic  coordinates  are  then  read  in  degrees,  and  changed  to  radians  by  a call 
to  CORRAD.  The  various  paths  from  transmitter  to  target  and  then  defined 
(RLA(ITCT),  RLO(ITCT)),  and  a call  to  subroutine  GEOPTS  defines  the  incremen- 
tal path  coordinates.  The  call  to  subroutine  GEODI  returns  the  total  distance  from 
transmitter  to  target  (SBKMS)  by  a Sodano  calculation.  This  value  is  later  used 
to  calculate  the  primary  wave  delay.  With  the  incremental  geographic  coordinates 
known  along  the  geodesic,  the  call  to  subroutine  GEORET  returns  the  ground  elec- 
trical properties  of  elevation  and  impedance  through  common  blocks /GROUND,  and, 
SDRDI/for  use  in  subroutine  INEQ.  Subroutine  INEQ  determines  the  secondary 
phase  factor,  TIMSAV.  The  time  of  arrival  (TOA)  for  a given  path  (ITCT)  is 
determined  from  the  relationship: 


TOA  (ITCT)  = ENC  * SBKMS  + TIMSAV 


(2) 


whore  ENC'  ground  refractive  index  divided  by  the  velocity  of  light.  ITCT  is 

tiien  incremented  for  a new  path  between  the  next  transmitter  in  the  chain  and  the  1 

same  receiver,  and  program  control  returns  to  statement  15  where  RLA(ITCT)  j 

and  fU.O(ITCT)  are  redefined.  The  above  process  is  then  repeated  until  three  ; 

TOA'S  are  calculated.  The  required  time  difference  (TDD  then  computed  from  the  ' 

relation:  ,* 

I 

t 

TDl  - EDI  + TC1A(2)  - TOA(l)  (3)  | 

where 

EDI  = emission  delary  for  slave  1.  [ 

TOA(2)  = Time  of  arrival  for  slave  1 at  receiver,  TOASl. 

roA(l)  = Time  of  arrival  for  Master  at  receiver,  TOAM.  j 

When  the  calculations  are  performed  for  airborne  locations,  the  secondary  phase  ' 

factor  contains  an  altitude  correction  derived  in  subroutine  INEQ  and  defined  as  ; 

ALTTMSV.  ; 

f 

Upon  completion  of  the  time  difference  calculation,  control  is  transferred 
back  to  statement  5 in  the  program  where  the  information  for  the  next  target  is 
read  in  and  the  entire  process  repeated.  The  program  exits  when  the  step  size 
(ADELS)  on  the  target  card  is  equal  to  or  less  than  zero. 

3.5.2  SUBROUTINE  CORRAD  I 

Subroutine  CORRAD  transforms  degrees  into  radians  for  a given  latitude  or  '• 

longitude.  The  subroutine  statement  is  SUBROUTINE  CORRAD  (RCOR,  IDEG, 

IMIN,  SEC,  ID,  IS,  lERR)  where: 

RCOR  = Location  in  radians. 

IDEG  = Location  in  degrees.  ' 

IMIN  = Location  in  minutes.  ' 

SEC  = Location  in  seconds. 

ID  = Character  for  north,  south,  east,  or  west.  [ 

IS  = Latitude  or  longitude  indicator.  I 

lERR  = Error  code.  i 

This  subroutine  is  called  from  the  driving  program  and  returns  radians  to  the  ! 

driving  program  through  the  argument  list.  It  is  used  to  transform  the  input  1 

transmitter  and  receiver  coordinates  into  radians.  | 
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3,5.3  SUBROUTINE  GEOPTS  j 

Given  the  distance  Sp  from  point  A to  a point  P on  the  geodesic  between  two  j 

prescribed  points  A,  B,  on  the  surface  of  a spheroidal  earth,  the  FORTRAN  sub-  j 

routines  GEOPTS  returns  the  latitude  0 , and  longitude  S of  P.  and  the  forward  I 

P P • I 

azimuth  (!p  of  the  geodesic  of  P as  shown  in  Figure  6.  j 


Figure  6.  Geodesic  Geometry 

The  derivation  of  the  equations  under  the  subroutine  are  described  by  Spies.  ^ 
The  subroutine  statement  is:  SUBROUTINE  GEOPTS  (SMP,  Rl^ATP,  RLONP, 
RAZP)  where 

SMP  - Sp  = distance  (in  meters)  of  P from  A, 

RLATP  = 0p  = latitude  (in  radians)  of  P, 

RLONP  * '^p  = longitude  (in  radians)  of  P, 

and 

RAZP  " " forward  azimuth  (in  radians)  of  geodesic  of  P. 

Latitudes  0^,  0g  and  longitudes  0^,  flg  of  the  prescribed  end-points  A,  B are 
stored  in  the  COMMON  block  PATH,  along  with  certain  output  parameters.  We 

6.  Spies,  K.  P.  (1975)  The  Analytical  Basis  of  Hufford's  Computer  Technique  for 
Determining  Topographm  Profiles,  Institute  for  Telecorn  munic  at  ion  Services 
Memo  dated  Nov.  4,  1975. 
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I 


I 
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have  the  statement:  COMMON/PATH/Pl.ATA,  HLONA,  HIATB,  RLONB.  HAZA, 
RAZB,  SMB,  where 

RLATA  - latitude  (in  radians)  of  A, 

RLONA  ~ - longitude  (in  radians)  of  A, 

RLATB  = 0p  = latitude  (in  radians)  of  B, 

RLONB  - fig  = longitude  (in  radians)  of  B, 

RAZA  ~ 't  p ~ forward  azimuth  (in  radians)  of  geodesic  of  A, 

RAZB  " = forward  azimuth  (in  radians)  of  geodesic  of  B, 

and 

SMB  ' ’ Istigth  (in  meters)  of  geodesic  from  A to  B. 

North  latitudes  and  east  longitudes  are  positive,  whereas  south  latitudes  and  west 
longitudes  are  negative.  Positive  azimuths  are  measured  clockwise  from  north. 

Parameters  which  specify  relevant  dimensions  of  the  spheroidal  earth  are 
defined  in  the  DATA  statement:  DATA  ASPH,  ESQ,  CESQ/6378206.  4,  0.006768658, 
0.993231342/  where 


ASPH  = a = length  (in  meters)  of  semimajor  axis  of  spheroid, 

2 2 2 2 

E.SQ  ■ e = (a  - b )/a  (where  b is  length  of  semiminor  axis  of  spheroid), 

CESQ  = 1 - e^, 

and  the  numerical  values  correspond  to  the  Clarke  spheroid  of  1866.  Any  other 
desired  spheroid  can  be  used  in  GEOPTS  simply  by  replacing  the  above  DATA 
statement  by  one  containing  the  appropriate  numerical  values.  For  example,  for 
the  International  spheroid,  one  would  use: 

DATA  ASPH,  ESQ,  CESQ/6378388. 0,  0.  006722670,  0.993277330/. 

As  we  shall  see  in  the  following  section,  the  quantities  0p,  5^,  il are  approx- 
imated by  esculatory  cubic  polynomial  functions  of  the  distance  Sp  of  P from  A;  for 
a fixed  spheroid,  the  coefficients  in  these  polynomials  depend  only  on  the  location 
of  the  end-points  A,  B.  To  achieve  computational  efficiency  in  situations  where 
0 , 9 , <l  are  to  be  calculated  for  more  than  one  point  P on  the  geodesic  between 
a fixed  pair  of  end  points,  GEOPTS  is  provided  with  a second  entry  point  PCOORD 
(by  including  the  ENTRY  statement  ENTRY  PCOORD),  whereby  the  subroutine 
skips  the  coefficient  calculations  and  proceeds  directly  to  the  evaluation  of  the  poly- 
nomial expressions  for  0p,  9^,  li For  a given  pair  of  end-points  A and  B,  the 
initial  call  to  the  subroutine  must  use  the  main  entry  point  GEOPTS:  that  is,  the 
calling  program  reference  must  be 
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L. 


CALL  GEOPTS  (etc.,  etc 


) 


Since  the  polynomial  coefficients  have  already  been  evaluated,  subsequent  calls 
(for  that  path)  should  then  use  the  entry  point  PCOOllD:  that  is,  the  calling  pro- 
gram reference  should  be: 

CALL  PCCKJRD  (etc.,  etc.,  . . . ). 

It  is  thus  evident  the  subroutine  GEOPTS  is  particularly  efficient  in  those  cases 
where  0p,  6^,  :/p  are  desired  for  several  to  many  points  P on  a single  geodesic. 

3.5.4  SUBROUTINE  GEORET 

This  subroutine  is  called  from  the  driving  program  LORANCO.  Subroutine 
GEORET  obtains  from  the  data  through  GETELV  the  values  of  elevation  and  imped- 
ance which  occur  at  intervals  along  a geodesic  specified  by  arrays  LAT  and  LON. 
The  path  elevation  and  impedance  data  are  transmitted  to  the  driving  program  via 
blocked  common  statements. 

The  subroutine  statement  is  SUBROUTINE  GEORET  (IJ)T,  LONG,  DSIvM, 

NP)  where: 

LAT  = Array  of  latitudes  along  a geodesic  path. 

LON  = Array  of  longitudes  along  a geodesic  path. 

DSKM  = Distance  along  geodesic  path  in  KJVl. 

NP  = Number  of  points  along  the  geodesic  path. 

To  obtain  elevations  and  impedances  from  the  data  base,  GEORET  steps 
through  NP  latitude  and  longitude  points  calling  GETELV  each  time.  Amplitude 
and  phase  data  are  returned,  converted  to  real  and  imaginary  values,  and  stored  in 
arrays  DR  and  DI.  Elevation  information  is  stored  in  array  Z.  Subroutine 
GEORET  then  calls  subroutine  INT  to  obtain  the  first  and  second  derivative  of  the 
elevation  data. 

Support  information  required  by  GEORET  and  supplied  by  common  blocks 
include: 

(1)  Common/GROUND/ 

(2)  Common/INDUCT/ 

(3)  Common/TE/ 

(4)  Common/SDRDI/ 

(5)  Common/CITCT/ 

By  use  of  common  blocks,  certain  of  the  information  can  be  directed  to  the 
double  integral  subroutine  (INEQ)  and  its  support  subroutines,  while  other  control 
data  are  transmitted  to  the  driving  section  of  the  overall  program. 


I 


j 3.5.5  SUBROUTINE  GETELV 

‘ This  subroutine  is  called  from  subroutine  GEORET.  Subroutine  GETEI.V 

P 

( indexes  the  input  latitude  and  longitude  coordinates  for  subroutine  UNPACK  to 

I obtain  the  proper  ground  electrical  properties  off  the  disc. 

: The  subroutine  statement  is  SUBROUTINE  GETELV  (LAT,  LONG,  AMP, 

FAZ,  ELVTN)  where: 

LAT,  LON  - The  input  coordinate  of  each  incremented  point  along  the 
geodesic  path. 


I 


I 


f 
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AMP,  FAZ,  ELVTN  = The  output  complex  impedance  (AMP,  FAZ)  and 
elevation  (ELVTN)  from  the  disc  data  base. 

Two  indices  are  selected  of  adjacent  longitude  records,  NPACK,  which  encom- 
pass the  input  coordinate.  These  two  records  are  separated  by  30  in.  in  longitude 
and  cover  a latitude  of  one  degree  or  120  points.  The  elevation  data,  ELVTN,  is 
obtained  by  a linear  interpolation  within  the  30  in.  square  surrounding  the  desired 
coordinate,  and  is  stored  in  eleven  bits  of  the  available  thirty  bit  word.  The 
impedance  data  of  the  southwest  corner  of  the  30-in.  square  is  used  to  represent 
the  entire  square.  No  impedance  interpolation  is  required.  For  the  complex 
impedance  data,  eight  bits  are  used  for  the  magnitude  and  eleven  bits  for  the 
phase. 

The  impedance  and  elevation  data  are  returned  to  subroutine  GEORET  through 
the  argument  list. 

3.5.6  SUBROUTINE  SETUP 

This  subroutine  provides  the  spheroidal  data  for  subroutine  GEODI.  Inputs 
are  the  semimajor  (AO)  and  semiminor  (BO)  axis  of  the  earth  in  meters.  The 
spheroidal  flattening  (FL)  and  eccentricity  square  (ESQ)  are  calculated.  The 


spheroidal  constants  for  various 

Ellipsoid 

ellipsoids  are: 

Semimajor 

Semiminor 

International 

6378388.  0 

6356911.  9 

Clarke  1866 

6378206.  4 

6356583.  8 

Clarke  1880 

6378249.  1 

6356514.  9 

Bessel 

6377397.  2 

6356079.  0 

Subroutine  SETUP  is  called  from  subroutine  GEODI  and  returns  the  required 
constants  through  the  COMMON/GENERA  L/block. 

3.  5.  7 SUBROUTINE  INT 

Subroutine  INT  is  called  by  subroutine  GEORET  and  calculates  the  first  and 
second  derivative  of  the  elevation  at  each  point  along  the  geodesic  path.  The  sub- 
routine statement  is  SUBROUTINE  INT  (I,  K)  where: 


1 
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I = Position  in  x and  z arrays  on  which  to  center  calculations. 

K = Position  in  arrays  Z and  ZP  to  store  calculated  values. 

An  analytical  expression  for  a parabolic  fit  to  three  data  points  closest  and 

7 

including  the  required  point  along  the  geodesic  path  is  derived  by  .Sheed.  Differ- 
entiating this  curve  twice  yields  the  required  first  and  second  derivative  of  the 
elevation  at  the  position  on  which  the  calculations  are  centered.  The  results  are 
returned  to  GEORET  for  use  in  subroutine  INEQ  via  Common/Ground/block. 

3.  5.  8 SUBROUTINE  UNPACK 

Subroutine  UNPACK  takes  each  sixty  word  record  from  the  data  base  and 
unpacks  each  word  into  two  complete  data  points  composed  of  elevation  and  imped- 
ance information.  Each  called  record  contains  the  data  for  120  points.  The 
original  sixty-bit  words  when  unpacked,  allow  thirty  bits  for  each  data  point. 

These  are  allocated  as  follows; 

Elevation  in  meters,  ELEV,  11  bits,  ±7  meters. 

Magnitude  of  Impedance,  AMP,  8 bits,  ±2  ohms. 

Phase  of  Impedance,  FAZ,  11  bits,  ±0.001  radians. 

The  information  obtained  is  returned  to  GETELV. 

For  machines  with  32  bit  words,  the  UNPACK  subroutine  is  altered  so  that  one 
data  point  is  obtained  from  each  word. 

3.  5.  9 SUBROUTINE  INEQ 

Subroutine  INEQ  is  called  from  the  driving  program,  LORANCO,  and  returns 
the  secondary  phase  factor  or  additional  time  correction.  The  LAT  and  LON  in 
the  calling  statement  are  the  array  of  path  length  latitudes,  and  longitudes 
respectively. 

The  following  constants  are  first  set  for  the  operation  of  this  program: 

J 

NUT  = 0 or  1 depending  on  type  of  receiver  antenna. 

RAD  = radius  of  earth  = 6.  36739  X 10®  meters. 

ALPHA  = Vertical  lapse  factor  = 0.  85. 

FREQ  = Frequency  = 0.  1 MHz. 

N.START  = Index  of  first  distance  at  which  the  field  is  to  be  found  as  a function 
of  altitude  = 0.  5 km. 

ETA  = Refractive  index  = 1.000338. 

7.  Sheed,  F.  (1962)  Theory  and  Problems  of  Numerical  Analysis,  Shaum's  Out- 
line Series,  McGraw  Hill  Co. 
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The  attenuation  function  is  calculated  at  discrete  points  as  a function  of  geo- 

O 

desic  distance.  The  amplitude  or  modulus  of  the  complex  attenuation  function  is 
the  field  intensity  and  the  phase  or  argument  is  the  secondary  phase  correction  in 
radians.  To  solve  for  this  attenuation  function,  the  boundary  conditions  must  be 
known,  then  the  solution  can  be  extended  step  by  step  by  numerical  integration. 

By  assuming  that  the  ground  is  smooth  just  in  front  of  the  transmitter,  the 

9 

first  few  points  can  be  calculated  with  classical  formulas.  For  the  remaining 
points,  an  integral  equation  approach  is  employed  which  allows  one  to  introduce 
variations  in  ground  elevation  and  variations  in  ground  impedance  relative  to  a 
classical  spheroid.  Special  care  must  be  exercised  in  the  imegration  near  the 
beginning  and  end  of  the  integration  path  because  the  integrand  approaches  infinity 
at  these  points.  A Gaussian  quadrature  integration  formula  is  employed  in  the 
area  of  such  singularities,  and  Simpson's  rule  is  used  in  the  rest  of  the  interval. 
The  effective  ground  impedance  which  combines  the  elevation  and  impedance  varia- 
tions is  obtained  through  a common  statement  from  subroutine  GROUND.  After 

the  calculations  along  the  ground  are  completed  for  a given  transmitter  to  receive 
4 11 

path,  a height  gain  loop  ■ is  activated  if  the  altitude  input  data  is  other  than 
zero.  Upon  completion  of  the  integration,  the  results  are  returned  to  the  driving 
program  through  common  block 'DELAY/. 

Input  information  initially  obtained  by  the  driving  program  from  the  stored 
disc  data  base  is  transmitted  to  subroutine  INEQ  via  common  blocks: 

Common/GROUND /and  Common/TE/ 

Transmitter-receiver  path  increment  data  enters  via  common  blocks: 
Common/SS/and  Common/SDRDI/ 

Integral  equation  control  of  variables  from  the  driving  program  into  subroutine 
INEQ  enter  via: 

Common/CITCT 

Calculations  of  secondary  phase  factors  in  units  of  microseconds  are  returned  to 
the  calling  program  via: 

Common/DELAY  / 

8.  Johler,  J.  R. , and  Horowitz,  S.  (1973)  Propagation  of  LORAN-C  Ground  and 

Ionospheric  Wave  Pulses.  Office  of  Telecommunications  Report  73-20 
(Superintendent  of  Documents,  U.  S.  Government  Printing  Office, 

Washington,  D.C.  20402). 

9.  Johler,  J.R.,  Keller,  W.J.,  and  Walters,  L.  C.  (1956)  Phase  of  the  Low 

Radiofrequency  Ground  Wave,  NBS  Circular  573  (Superintendent  of  Docu- 
ments, U.  S.  Government  Printing  Office,  Washington,  D.C.  20402 

10.  Abramovitz,  M.  and  Stegun,  I.  (1964)  Handbook  of  Mathematical  Inunctions  — 

NBS  — Applied  IVlath  .Series  55,  U.  S.  Government  Printing  Office. 

11.  Scott,  R.  (1966)  Phase  of  the  Height -(Sain  Function  of  the  Lov.  I'requency 

Ground  Wave,  Report  2900-156-T  of  Project  Michigan  (Willow  Run 
Laboratories,  Ann  Arbor,  Michigan  48106). 
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II. I)  I’riifimni  for  I’iirl  III 


PkOGF  if-  LCF/lNCn(  lNPUT  = l£e  , CUTP0T  = 1^  8 ,T  APE  2) 

Dlf-FUSION  irENTF(3) 

DIKfNSTOF  lAT (999) ,L0N(9g9) 
iJlMENSION  CIETUUMIl)  ,TIMC'U*<(3) 

DIMENSION  TOA(3) 

OIPENSTON  CTSTSOO(3) 

OIMENSICN  TFH(3) 

GlhENSIOH  ALITt.S\/<3,2) 

DIMENSION  FAIHL(3) 

DIMENSIC  f LCUM  (3  ,t)  , AOOM(  3,<.),RLA(3),'»LO(3) 

CDMMCN/F  ATF/RLAT A, RLONA, PLATB, FLONP , PA ZA i RA ZB , SBM 
COMMCF /CITCT/IICT,LTOP,ALTT  MSV 
CDMMCN/CEL  AY/NP, DISTSA V,T IMEA V 
CO MM CN/Z OTA ZAP PAY (AD) 

NAMFLIST/NAMl/7  0A,TCil,TD2  iFrjI,E02,0IST0UM,TIMnUM,0ISTSC0 
1 , T P W 

DAI  A FTODEG/57.  2997795130S23Z 
DATA  ITCT/1/ 

DATA  APRAY/lO0.,1.300338,‘*C.,.l,5.,l.,3.,3G00.,lC.,C.,lOC0., 

1 •89tC«fl«f2  6*  C • / 

C=. 299792 

EN:l.tt,C338 

ENC=FN/C 

C RcAO  TRANSMITTER  COOROINAT- S FOR  MASTER,  S1,S2. 

PRINT  101 
PRINT  13 

101  FORMAT  (♦ITPANSMITTFR  COORDINATES*) 

MIbKOO=C 

DO  913  1=1,3 

READ  1C2  ,IDEMTF,LOUM  (1,1)  ,L  CUM  (1,2)  , AOUMd,  1)  , ADUM  (1 , 2)  , 

1 LOUM(I,3) ,LbUM(I,M),ACUM(I,3),A0UM(I,A),TCA(I) 

PRINT  1C  3,ICENTF  ,LDUM(I,1)  ,LOUH  (I  ,2)  ,AOUM  (1,1), AOUMd, 2), 

1 LDUM(I,3),LDUM(I,A),ADUM(I,3) , ADUM (I, A) 

102  FOf.  MAT(2X,3A8,2(I5,I3,F7.3,fll)  ,F12.  1,FE.2,I3) 

103  FOPMAT(2X,3A8,2(2X,I9, I3,F7.3,Al) ) 
im  FOKMAT(*lhECEH/LR  COOKOINATFS*  ) 

CALL  CCPPAC(PLC(I),LDUM(I,3) , LOOM (I, A) .ADUM (I, 3) , A DUM ( I , A ) , 1 , M ISKO 
ID) 

call  CCRPACfPLA ( I)  ,LDUM(I,l),LDOM(T,2) , ADUM (1,1) , A CUM (1,2) ,0,MISKO 
ID) 

IF (MISKOC.  EO. 0 ) GO  TO  913 
PRINT  23 
GO  TO  10 
913  CONTINUE 
ED1=TCA (2) 

C02=TCA(3) 

5 CONTINUE 

READ  102,  ICEt'TP,LATDEGB,LATMINB,SECLATP,L  ATI0  6,LONC2  GE  ,LONMINe  , 
ISECLONE,!.  0MDB,AlT  ,ADELS,  NUT 
IF(ADELS.LE.L)  GO  TO  10 
PRINT  IOh 
PRINT  13 

PRINT  It  3, IPENTF ,L ATO£Ge,LATMINP,SECLATB,LATiDB,LONDFGE,LCNMINe, 
1SECL0NB,LCMDB 
A?r  AY  (10  ) = ALT 
ARRAY  (lA ) = NUT 
MISKOO  = C 


I 

► 
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CALL  C0RRAG(RLATa,LATDEGB.LATHIMB,SECLAT9.LATIoe.0,HISK0D) 

CALL  CORO AOIRL 0MB, LONDE&B,LONHIN8,SECLONe,LtmiDB>lt>1ISKOD) 
IF(NtSKOO.EO.O)  60  TO  15 

29  PitINT  23 

23  rORMAT(//llX,  3H***,  9X,  29HERROR  IN  EMO-POIMT  COORDINATE,  lOX,  3H 
♦•••/IlX.SAH***  CALCULATIONS  FOR  THIS  PATH  HAVE  BEEN  DELETED  •••) 
60  TO  5 
15  CONTINUE 

RLATA^RLAdTCT) 

RLONA>RLO(ITCT) 

13  FORMAT (12H0  LOCATION,22X,8HLATITUDE,9X,9HLONGITUOE/ 

1 26X,2(5X,13HI0EG-HIN-SEC)>/) 

LTOP»0 

25  CALL  GEOPTS(0.0,RLATP,RLONP,RAZP) 

CALL  6EOOI(RLATA,RLONA,RAZAS,SeHS,RLATe,RLONB,RAZeS) 

S3KNSsSBNS«.e01 

OAZA  = RT00E6«RAZA 

OAZB  = RT00E6*RAZB 

SBKH  = SeM*1.0E-3 

OELaSBKN'SeKNS 

PRINT  2, ITCT,SBKHS,SBKM,OEL 

2 FORMAT t*OGEOOESIC  DISTANCES,  TRANSHITTER* I 2*  IPOINT  A)  TO  RECEIVER 
1 (POINT  e)*/7X*SOOANO*9X»HOFFORO*6X*DIFFERENCE*/3F15.5) 

NP  s SBKH/AOELS*1.0 

OSKN  = NP 

OSM  s SBN/DSKH 

OSKM  = DSM*1.0E-3 

PRINT  33  , SBKH, DSKN, OAZA, DAZ6 

33  FORHATI/ZiaX,  20HLENGTH  OF  GEODESIC  =,  F11.5,  IIH  KILOHETERS/ZieX, 

♦ ZSHOISTANCE  INCREMENT  -,  F11.5,  IIH  KILOMETERS/// 12X , 3AHF0RMA*0 
♦AZIMUTH  OF  GEODESIC  AT  A s,  F12.6,  BH  DEGREES //12X , 3AHF0RHARD  AZI 
♦NOTH  OF  GEODESIC  AT  B *,  F12.6,  8H  DEGREES) 

LIT (1><L0UH(ITCT,1) *10600 ♦L0UM(ITCT,2)  *10 0*IFIX(  AOUM(  ITCT ,1)  ) 

LON(1)=LOUH(ITCT,3)*10000*LDUH(ITCT,4)*100»IFIX(AD UM(ITCT,3) ) 

NPP2=NP«-2 

SPH  - 8.0 

DO  100  IP=1,NPP2 

SPH  * SPN+OSM 

CALL  PCOORO(SPH,RLATP,RLONP,RAZP) 

HISKOD  = 0 

CALL  CORDHS(RLArP,LATOEGP,LATNINP,SECLATP,LATIOP,0 ,MISVOD) 

CALL  CORONS(RLONP,LONOE6P,LONHINP,SECLONP,LON1DP,1, HISKOD) 

IF  (HISKOD)  30,35,30 

30  PRINT  <*3,IP 

A3  F0RMAT(//2X,  62H***  ERROR  DETECTED  DURING  CONVERSION  OF  COOROINAT 
♦ES  FOR  POINT,  15,  5H  ••*/2X,  2H**,  TX,  52HFR0M  RADIANS  TO  (ALPHAM 

♦ ERIC)  OEGREES-MIHUTES-SECONDS,  7X,  3H***/2X,  3H*** , 66X,  3H***/2X, 

♦ 3H***,  7X,  52HFURTHER  CALCULATIONS  FOR  THIS  PATH  HAVE  BEEN  DELETE 
♦0,  7X,  3H***) 

GO  TO  5 
35  CONTINUE 

LATdP^l  )=LATOEGP*ll»00^LATNINP*100^IFIX(SECLATP) 

LONdP^l  )^ONOEGP*10000^LONNINP*100^TFIX(SECL0NP) 

100  CONTINUE 
IP=NP^3 

CALL  GEORET(LAT,LON,OSKN,IP  ) 

C ADJUST  NO  OF  PTS  TO  HATCH  INEQ2E. 
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NPsNP+1 

CALL  INEOZECLATtLON) 

T1)A(ITCT)  = ENC*S0KNS  ♦TIHSAV 

TPH(1TCT)=ENC*SBKMS 

DISTSOOIlTCT)»SBKI1S 

PATHL (ITCT)=SBKHS 

DIST0UIHITCT)=01STSAW 

TIMOUH(ITC1)=TIHSAW 

irCT=ITCT*l 

IF(ITCT.LT.i»)GO  TO  15 tlF I LTOP. LT.  1)  PRINT  AlO 
AlO  F0RHAT(*0  GROUNOHAVE  TINE  DELAY  CALCULATION*) 
T01*ED1*TOA(E)-TOA(1) 

T0Z=ED2*T0A(3)-T0A(1) 

PRINT  NANI 

IF(LTOP.LT.l)  GO  TO  6 
00  8 L=1,LT0P 
PRINT  7 

7 FORMAT(*0  ALTITUDE  TINE  DELAY  CALCULATIONS*) 

DO  9 1=1,3 

T0A(1)=ENC*PATHL(I)*ALTTHSV(I,L) 

9  CONTINUE 

T01=EOl*TOA(Z)-TOA(l) 

TD2=ED2+T0A<3)-T0A<1) 

PRINT  NANI 

8 CONTINUE 
6 CONTINUE 

ITCT=1 
GO  TO  5 

10  CONTINUE 
PRINT  11 

11  FORMAT(*C  END  OF  INPUT  DATA.*) 

END 


SUBROUTINE  GEOREKL  AT,LON  ,OSKM,NP) 

OIHENSION  LAT(999) ,LON(999) 

COHHON  /GROONO/ZOUH,  7,  ZPOUH,  ZP,  XDUH,  X,  RO,  DEN,  R,  U,  ONX,  MA 
IVE,  R2,  ZPP,  0,  DOOl,  0002 
CONHON/INOLCT/HAV 

COHHON/TE/TLA,TLON,RLAT,RLON,HLAT,HLON,  ZPPl, 0INC1,TPI, NPTS 

1,ZNAZ,HTH 

COHHON  /SDRDI7S(999) , 0R(999),  01(999),  LLH 

DIHENSION  20UN(6),  Z(1009),  ZPDUH(6),  ZP(1009),  XDUH(6),  X(1089), 

1 ZPP1(1009) 

COHHON/CITCT/  ITCT 
CONPLEX  C,CEN, 0001, 0002 
TPI=6,2e  3185307 
TlA=LAT(1) 

TL0N=L0N(1) 

RLAT=LAT (NF) 

RLON=LCN(NF) 

DINC1=0SXM 

NPTS=NP 

S(1)=0. 

X(1)=0. 

DO  1 1=1, NP 

AAHP=.07AAJFAZ  = .7A10  9$Ztl  )=2. 


CALL  GETELM(LAT(I) «LON(I) fAANP,FAZ,Z(I)  ) 

IF<I.6T.l)  :>tI)*StI-l>»DSt(N 
Oll(I)  > AAHP  « COS(FAZ) 

OKI)  = AAPP  • SIN(FAZ) 

CONTINUE 

LlH=NP 

GENERATE  INDEPENDENT  PATH  TERRAIN  ELEVATIONS  RELATIVE  TO  TRANSMIT 
QcZCl) 

Z(1I>:0. 

DO  2 1=2, NP 

FOR  NOW,  CONVERT  DISTANCE  TO  METERS. 

S(I)=S(I)«1000. 

X(I)=S(I) 

Z(I)=Z(I)-0 

CONTINUE 

OBRATN  ELEVATION  OERAVITTVES. 

NPN1*NP-1 
DO  3 L=2,NPM1 
CALL  INTCL.L) 

ZPPl  (L)*ZPP 

CONTINUE 

ZPfl)*0. 

ZPP1(1)*0. 

ZP(NP)=ZF(NP-1) 

ZPPl <NP)=ZPP1(MP-1> 

RETURN 

END 


SUBROUTINE  CORRAD(RCOR,IDEG, IHIN,SEC,IO, IS.IERR) 

DIMENSION  IDS (4) 

DATA  riDS=lHN,lHE,  IHS.IHH) 

ISS  = IS 
IF  (ISS)  10,5,15 
5 IF  (ID.EQ.IOS(l))  GO  TO  25 

IF  (IO.EQ.IDS(3))  GO  TO  30 

10  lERR  = 1 
RETURN 

15  IF  (ISS-1)  20,20,10 

20  IF  (IO.EO.IOS(2))  GO  TO  25 

IF  (I0.E0.IDS(4))  60  TO  30 

lERR  » 1 
RETURN 

25  SIGN  = 1.0 
GO  TO  35 
30  SIGN  = -1.0 
35  IF  (IOEG-160)  40,40,10 
40  IF  <IHIN-60)  45,10,10 
45  IF  (SEC-60.0)  50,10,10 

50  RCOR  = SIGA*FLOAT(IOEG)*(1.74532S25199433E-2)*FLOAT(IHTN)* 
> (2. 9 OSae 2006657 22E-4)»SEC«4. 048 136611 09536E-6 
RETURN 
END 
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SUBROUTINE  6E0PTS(SNP,RL*TP,RL0MI>.RAZP) 
CaHNON/PATH/RLATAtRLONA,RLATBtRLONB,RAZ«,RAZ8,SNB 
DATA  ASPHtESOt CESQ/63783eB.00l.  6«72267e022E>3,  S. 932773300E-1/ 
DATA  PI02f  Ply  TNOPI/t.57079632679i»90y  3. 1A15 926535«97«t 
* S.2831«53a717959/ 

IF  (ABS< PLATA) -PI02)  SylOylO 
5 IF  IABS(RLATB)-PT02)  ISylOylt 
10  PRINT  SyRLATAyRLATB 

3 FORNATIlHl,  69H  SUBROUTINE  GEOPTS  CALLED  WITH  ENO'POINT  LATITUO 
»E  OUT  OF  ACCEPTABLE  RANGE  <-PI/2y PI 72) //1 3Xy  16HLATITUOE  OF  A =y 
»FS.2y  29N  RADIANS  LATITUDE  OF  B >,  F0.2(  SH  RADIANS) 

PRINT  13 

13  FaRNAT(33X,  28HPROGRAH  EXECUTION  TERMINATED) 

STOP 

15  IF  (A6S(RL0NA)-TH0PI)  20y2Cf25 
20  IF  (AeS(RLONB)-TMOPI)  38y30y25 
25  PRINT  23,RLONA,RLONS 

23  FORHATdHly  OBH  SUBROUTINE  GEOPTS  CALLED  MITH  END-POINT  LONGITU 
«OE  OUT  OF  ACCEPTABLE  RANGE  (-2PI y2PI)//10Xy  17HL0NGITU0E  OF  A =y 
♦F9.2y  31H  RADIANS  LONGITUDE  OF  B xy  F9.2»  8H  RADIANS) 

PRINT  13 
STOP 

30  ALONR  < RLONA 
BLONR  s RLONB 
IF  (ALONR^PI)  35«35.90 
35  ALONR  * ALONR«TMOPI 
90  IF  (BLONR«PT)  9Sy95y50 
95  BLONR  X BLOMt^TMOPI 
50  BLON  X 8LONR-ALONR 
IF  IBL0N4PI)  55y70y75 
55  IF  (BLONR)  60y65y65 
60  BLONR  X BLONR»THOPI 
GO  TO  95 

65  ALONR  X ALOMt-TMOPI 
GO  TO  95 
71  PRINT  33 

33  FORMAT (IHly  88H  GEODESIC  PATH  INCLUDES  A GEOGRAPHIC  POLE  - SUBR 
♦OUTINE  GEOPTS  CANNOT  HANDLE  THIS  CASE) 

PRINT  13 
STOP 

75  IF  (BLON-PI)  95y70y80 
80  IF  (BLONR-PI)  90y9ly85 
85  BLONR  X 8LCNR-TMOPI 
GO  TO  95 

90  ALONR  X ALONR»TNOPI 
95  SMB  X RLATAHRLATB 
HETA  X SIN(0.5»SNB) 

HETA  X ( l.•-€SO•HETA•META)7CESQ 
HLAT  X ASP)‘/(HETA*S()RT(HETA*CESQ)) 

BLON  X BLONR-ALONR 
ALON  X HETA«8L0N 
Q X SIN(ALON) 

BLAT  - COSIRLATB) 

ALAT  X COS(RLATA) 

SAZA  X Q*BLAT 
SAZB  X QxalAT 
0 X SIN(0.5«ALON) 
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105 

110 

115 

120 


Q = Q*Q 

AAZ  « RLATB-RLATA 
CAZB  = (1.0-Q»*StN(AAZ) 

Q = Q*SIN(SHB) 

CAZA  = CAZE«Q 
CAZB  * CAZe-Q 

Q s SORT(SAZA«SAZA«CAZA»CAZA) 

SMB  s HLAT«ASIN(Q) 

SAZA  = SAZA/Q 
CAZA  X CAZA/Q 
SAZB  X SAZe/Q 
CAZB  X CAZB/Q 
RAZA  X ATAN2(SAZAtCAZA) 

RAZB  X ATANZtS AZB«CAZBI 

ALON  X RAZE-RAZA 

IF  ( AeS(ALCN)-PI)  120.105.105 

IF  (ALCN)  110.110.115 

RAZB  X RAZB^TMOPI 

GO  TO  120 

RAZA  X RAZA'fTMOPI 

H.AT  X CESC*ESQ*ALAT*AUAT 

HLON  X ASPH/SQRT(HLAT) 

HLAT  X HtOR*CESQ/HLAT 

HLON  X HLON* AL AT 

CLAT  X CAZA/HLAT 

CLON  X SAZA/HLON 

HLAT  X CES0*ESQ*8LAT*BLAT 

HLON  X ASPH/SQRT CHLAT) 

HLAT  X HLON*CESO/HLAT 
HLON  X hLON*BLAT 
ALAT  X CAZE/HLAT 

8LAT  X <3.0*AAZ/SHB-ALAT-2,0*CLAT)/SMB 
ALON  X SAZe/HLON 

BLON  X (3.0*BLON/SHB-ALON>2.0*CLON)/SHB 
AAZ  X 3.0*SHB 

ALAT  X ( (ALAT-CLAT)/SNB-2.0*BLAT)/AAZ 
ALON  X ( (ALON-CLON) /SNB-2.0*BLON) /AAZ 
AAZ  X (RAZB-RAZA)/SNB 
ENTRY  PCOORO 

RLATP  X ( (ALAT*SNP*BLAT)*SHP*CLAT)«SHP«RLATA 

RLONP  X ((ALON*SHP*BLON)*SHP>CLON)*SHP*ALONR 

RAZP  X AAZ*SNP«RAZA 

RETURN 

ENO 


SUBROUTINE  COROHS (RCOR.IDEG. ININ .EEC.IO .IS.IERR) 
DIMENSION  IOS(%) 

DATA  (IOSxim.lHE.lHS.lHH) 

RANG  X RCOR 

SEC  X AeS(RANG>*20  626<».a0  62<»70  96 
ISS  X IS 


IF 

(ISS)  5.10.15 

5 

lERR  X 1 

RETURN 

10 

IF 

(SEC-32A000.0005) 

25,5,5 

15 

IF 

(ISS-1)  20,20.5 

20 

IF 

(SEC- EAOOOO. 0005) 

25.5.5 

L. 


25  IF  IRANG)  30»35t35 
30  ISI  « 2 
GO  TO  <>0 
35  ISI  « 0 
40  lOEG  3 SEC/3600. 0 

ININ  > SEC/60. 0-60. 0*FLO«T(IOEG) 

SEC  » SEC-3600. 0»FLOAT(IDEG)-60.0*FLOATCIHIN) 
ISEC  ® SEC 

SEC  = SEC-FLOAT (ISEC) 

IF  (SEC-0. 9995)  60,45,45 
45  SEC  = 0.0 
ISEC  = ISECn 
IF  (ISEC-6Q)  60f50t50 
50  ISEC  « 0 

ININ  = IHINi-l 
IF  (ININ-60)  60.55.55 
55  ININ  = 0 

IDEG  = lOEGH 
60  LDX  = ISI^ISS^l 

BEC=FLOAT(ISEC) ♦SEC 
IO=IOS(LOX) 

RETURN 

END 


SUBROUTINE  INEQ2E( LAT .LON ) 

DIHENSION  LAT(999)fL0N(999) 

COHNON  /S5/0ISS(999  ) 

CONHON  /GROUNO/ZOUN.  Z«  ZPOUN.  ZP.  XOUN.  X.  RO.  DEN,  R.  U,  ONX.  NA 
IVE.  R2.  ZRP.  D«  0001,  0002 
COHHON/IhOtCT/HAtf 

CONNON/TE/TLA,TLON,RLAT,RLON.HLAT,HLON.  ZPPl.OINCl.TPI, NPTS 

l.ZHAZ.HTH 
COHNON  /ZOTA/ARRAV 

COHNON  /SDRDI/S(999) . 0R(999),  01(999),  LLN 
COMHON/CITCT/ITCT,LTOP,ALTT  HSW 
COHNON/OELAY/NP.OISTSAV.TINSAV 
OINENSION  ALTTNSV(3.2) 

DIMENSION  E(999  . 2).  F(999  , 2) 

OINENSION  ZZ(5).  H(5) 

OINENSION  Z0UH(6).  Z(1009).  ZPDUH(6>,  ZP(ie09).  X0UH(6).  X(1009), 
ITZER(l).  T(IOOO).  GM(5).  GX(5) . EM(9),  EX(9),  NG(5),  XS(9).  ME(9), 
2 N(IOOO) ,ARRAV(4t) ,ZPP1(19B9) 

COMPLEX  0,  OEN,  0001,  0002,  FLEAF,  FI,  F2,  F3.  Gl,  G2,  G3,  SUN,  TE 
IRN,  M,  ME.  MG.  MM.  MX 
COMPLEX  TS.FINO.FINOH 

EQUIVALENCE  (ARRAV(l).  FKHZ) . (ARRAV(2),  ETA),  (ARRAV(3),  DNAX),  ( 
lAltRAY(4),  OINC),  (ARRAY(8),  ZMIN)  , (ARRAY(9),  ZINC),  (ARRAY(IO).  Z 
2HAX).  (ARRAY(ll),  FLAT) 

DATA  lARRA/O/ 

DATA  ((GX(K).  K s 1,  5)  ■ .02216356001,  .1070315676,  .4615973615, 
1.7403346204,  .9404939262) 

DATA  ((GM(K),  K > 1,  5)  « .5910404494.  .5305334306.  .430172725. 

1 .2909026903,  .1333426006) 

DATA  ((EX(K),  K s 1,  4)  * .1051402026,  .3762245145.  .6909400124, 

1 .9  373342493) 
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DATA  ((EMCK),  K > !•  4)  .06568«519e9.  .196096265S,  .2525273456, 

1.1523625357) 

DATA  (TZER(1)«1.0) 

6  FORMAT  (//E20.7,5X,«ARRAVI1>  = F»«Z  = FREQUENCY  IN  KHZ*/ 

1 E20.7,5X,«ARRAr(2)  > ETA  s air  INDEX  OF  REFRACTION  AT  GROUND*/ 

2 E20.7,5X,*ARRAV (3)  = DNAX  = MAXIMUM  DISTANCE  OETHEEN  TRANSMITTER* 

3 ••  AND  RECEIVER  IN  KM*/ 

4 E20.7,5X,*ARRAV(4)  = DINC  = DISTANCE  INCREMENT  IN  KM*/ 

5 E2l.7,5Xy*ARRAV(5)  = NSTART  » INDEX  OF  THE  FIRST  DISTANCE  AT*, 

S • NHICH  THE  FIELD  IS  TO  BE  FOUND  AS  A FUNCTION  OF  ALTITUDE*/ 

7 E2S.7,5X,*ARRAV(6)  = NZINC  = THE  FIELD  IS  TO  BE  FOUND  AS  A *, 

8 •FUNCTION  OF  ALTITUDE  EVERY  NZINC  INCREMENTS  IN  THE  •/ 

9 40X,*DISTANCE  ARRAY  ••  NZINC  MUST  NEVER  BE  LESS  THAN  1*1 

7 FORMAT  (E20.7,5X,*ARRAV(7)  = INDEX  OF  THE  LAST  DISTANCE  AT  WHICH*, 

1 * THE  FIELD  IS  FOUND  AS  A FUNCTION  OF  ALTITUDE*/ 

2 E20.7,5X,  *ARRAT(B)  * ZMIN  = THE  MINIMUM  ALTITUDE  ABOVE  THE  *, 

3 'SURFACE  AT  WHICH  THE  FIELD  IS  FOUND,  METERS*/ 

4 E20.7,5X,  'ARRAY (9)  = ZINC  = THE  ALTITUDE  INCREMENT,  METERS*/ 

5 E20.7,5X,*ARRAY(10)  = ZMAX  s THE  MAXIMUM  ALTITUDE,  METERS*/ 

6 E20.7,5X, 'ARRAY (11)  = FLAT  = THE  FACTOR  USED  IN  THE  FLAT  EARTH  *, 

7 'TEST,  METERS*/ 

8 E20.7,EX,*ARRAY(12)  = ALPHA  = VERTICAL  LAPSE  FACTOR*/ 

9 £20. 7, 5X, 'ARRAY (13)  = EFFECTIVE  EARTHS  RADIUS*/ 

A E20.7,5X,*ARRAY(14)  = NUT*//) 

IF(IARRA.GT.l)  GO  TO  9999 
IARRAX2 

C SET  CONSTANTS 
NUTsARRAY(14) 

RAO  X 6. 3e  739E*6  / ARRAYdZ) 

ARRAY(13)  X RAO 
FREO  * FKHZ  * l.E-3 
NSTART  = ABRAY(5) 

NZINC  X ARRAV(6) 

IF  (ARRAV(6)  .LT.  1.)  NZINC  x i 

NZEROxO 

T(NZERO)xi,0 

WAVE  X 2,0'5e44729E-2  * FREQ 
WAVE  X WAVE  * ETA 
NAVxWAVE 

AMICRO  X 1,0  / (TPI  • FREQ) 

TX  X SORT(FREQ  * ETA)  * .0408389549 
X(l)  X 0. 

W(l)  X 1. 

9999  CONTINUE 

DECTX»X(2) 

DMAXxX(NPTS)/1000. 

DINCxOECTX/1000. 

ARRAY(4)=OINC 

OINClxOINC 

ARRAY  (7)  X ohAX  / DINC  ♦ 1. 

NZENO  X array (7) 

PRINT  6,  (AI:RAY(I),Ix1,6) 

PRINT  7, (ARRAY(I) ,1x7,14) 

NOFLAT  X 1 
NGO  X 1 

0 X CMPLX(OR(l),  OKI))  * SORKFREQ  / 0.1) 

INOST  X DMAX  / DINC  * .01 
DELTX  X 1000.  • DINC 
T(l)  X 1.  / SORT (DELTX) 
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I 


1 

I 


31  DO  32  K = 1,  5 

X(  - K»  = OELTX  • GX(K) 

CALL  GROUND  ( -K,  2,  0,  0) 

32  NG(K)  3 FLEAFIHA^Ef  O.t  O.t  X(  - K) , D) 

1 = 2 

37  CONTINUE 

ZPPsZPPl <11 

TCI)  = 1.  / SORTCXd)  ♦ DELTX) 

OMX  = OMCOS<X(I)  7 RAD) 

IGO  = 2 - HOOCI,  2) 

IL  = I - ICO  - 1 

RO  = SORT  (2.  • RAD  • (RAD  ♦ Z(I))  • OHX  ♦ ZCI)  • • 2) 

GO  TO  (40,  45) I NOFLAT 

40  IF  (I  .LE.  4 .OR.  (FLAT  .GT.  X(I)))  41,  45 

41  H(I)  = FLEAF(HAVE,  0.,  0.,  X(I),  D) 

GO  TO  90 

45  SJH  = 0. 

NOFLAT  = 2 
4S  00  50  K = 1,  5 

CALL  GEOHd,  • K,  1,  0.) 

TEKH  = HG(K)  » C«PLX(COSt NAVE  ♦ R),  - SIN(NAVE  * R) ) • DEN 
50  SUH  * (U  • GN(K))  7 SaRT(X(I)  - X(  - K))  • TERN  ♦ SUM 
SJM  = 3.  • T(l)  * SUM 
KK  = 1 

IF  IIL  .LT.  3)  GO  TO  100 

DO  60  K = 3,  IL 

CALL  GEOMd,  K,  1,  0.) 

TERM  = U * T(K  - 1)  • N(K)  • CMPLX(COS (NAVE  • R) , SIN(  - NAVE  • R) 
1)  • DEN 

GO  TO  (53,  55),  KK 


53 

SJM  = 4. 
K<  = 2 

60  TO  60 

• T(I  - K) 

• TERM  ♦ 

SUM 

55 

SUH  X 2. 
KX  X 1 

* Td  - to 

• TERM  ♦ 

SUH 

60 

CONTINUE 

100 

CONTINUE 

CALL  GEOMd,  2,  1,  0.) 

SUM  = T(I  - 2)  • T(l)  * U • M(2)  • CMPLX(COS(MAWE  • R),  - STN(NAVE 
1 • R) ) ♦ DEN  ♦ SUM 

62  CALL  GEOMd,  I - IGO,  1,  0.) 

F2  = U * T(IL)  • Nd  - IGO)  • CMPLX  ( COS  ( NAVE  • R),  - SIN(NAVE  • R) 
1)  • DEN 

SUM  = (SUM  ♦ TdGO)  • F2)  • .333333333  • OELTX 

GO  TO  (65,  66) , IGO 

65  FI  » F2 

F2  = TERM 

GO  TO  70 

66  CALL  GEOMd,  I - 1,  1,  0.) 

FI  = U * T(I  - 2)  • Nd  - 1)  * CMPLX  (COS  (NAVE  • R)  , - SIN(MAVE  • R 
1))  • DEN 

SUN  = SUM  ♦ .08333333333  • OELTX  • (5.  • T(l)  • FI  ♦ 8.  • T(2)  • F 
12  - T(3)  • TERM) 

70  Q * TX  7 Td  - 1) 

RHO  = 1.  ♦ 7P(I)  • • 2 

RNO  = ZPP  7 RHO 

75  TERM  r 1.2  7 T(l)  • Td  - 1)  • CMPLX(Q,  Q)  • (0  ♦ CMPLX(0.,  - RHO 
17  NAVE)) 

NX  = 1.  - CMPLX(Q,  0)  • (SUM  ♦ 2.  7 Td)  • (.466666667  • FI  - .066 
1666667  • F2))  7 (1.  ♦ TERM) 
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S5  M(l)  = HX 
90  OIST  = .001  * X(II 

CALL  INOF(0.,X(I) >NUT,FINO) 

TS=HCI)*FIHD 

AMP=CA8S(TS) 

PHArCANGCTS) 

PHIC  = - (PHA  - HAVE  • (RO  - X(T))) 

SEC  = PHIC  ♦ AHICRO 
IF(I.Nt.NP)  SO  TO  9076 
DISTSAV*OIST 
TIMSAV=SEC 
9876  CONTINUE 
I = 1 ♦ 1 

IF(I.LE.NPTS)  60  TO  37 
150  NGO  = Z 

IF  (7NAX.LE.0.)  RETURN 

irOP  = (NZENO  - NSTART)  t N7INC  ♦ 1 

IF  (ITCP  .GT.  1500)  GO  TO  280 

LT0P« (ZMAX-ZMIN) 7ZINC+1.1 

IF(LTOP.CT.5)  GO  TO  280 

FZ=(0.,0.) 

G2=(0..0.) 

OLDR^O. 

IF  (ITCP  .IT.  1)  GO  TO  105 
NPALT=NP-NSTART*1 
03  200  IS  = 1,  ITOP 

I = NSTART  ♦ (IS  - 1)  • MZINC 

II  = I 

IF  (I. GE. NZENO)  II^NZENO-l 

ISO  = 2 - POOd,  2) 

00  155  K = 1,  4 

XdOOO  ♦ K)  = X(I)  - OELTX  • EX(K) 

XS(K)  = XdOOO  * K) 

155  CALL  GROINO  ( 1 OOO^^Kt It.O,  0) 

CALL  CNEVKENCX,  H,  INOST,  XS,  WE.  4.  5.  1) 

IF  (LTOP  .LT.  1)  GO  TO  110 
DO  195  L = 1,  LTOP 
H(L)  = ZMIN  ♦ (L  - 1)  • ZINC 
IF  (H(L)  .LT.  ZCD)  156.  154 

156  WW  = 0. 

GO  TO  193 

154  ZZ(L)  = RAO  ♦ H(L) 

RO  = SQRT(2.  • RAD  • ZZIL)  • ONCOS(X(I)  / RAO)  ♦ H(L)  • • 2) 
SJH  = 0. 

00  157  K = 1,  5 

CALL  GEOMd,  - K,  2,  H(L)) 

TERM  = WG(K)  * CMPLXt COS (WAVE  • R>.  - SIN(MAVE  • R))  • OEM 

157  SUM  = U • GW(K)  • TERN  / SQRT(X(I)  - X(  - K))  ♦ SUM 
SUM  = 3.  *1(1)  • SUN 

CALL  GEOMd,  3,  2,  H(L)> 

FI  = T(l)  • Td  - 2)  * U • DEN  • W(2) 

G1  = CMPLX(COS(WAVE  • R) , - SIN(WAVE  • R) ) 

SUM  = SUM  ♦ FI  • G1 
ILO  = I - IGO 
KSO  E 1 

IF  (ILC  .LT.  3)  GO  TO  115 
00  182  K 3 3.  ILO 
F3  = F2 
63  = G2 
F2  = FI 
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I 


GZ  > G1 

CALL  GEOHCI,  K,  Zy  H(L)) 

G1  » CHPLX  (COS (HAVE  • R) , - SIN(MAVE  • R) J 
FI  = T<K  - 1)  • T(I  - K)  • U • OEM  • M(K) 

OELTG  = HAVE  • (R  - OLOR) 

GO  TO  (158»  159.  172)  . KGO 
15*  SDH  « SOH  ♦ 4.  • FI  * G1 
KGO  s 2 
GO  TO  102 

159  SJH  = SUM  ♦ 2.  • FI  • G1 
KGO  s 1 

IF  (AeS(OELTG)  .GT.  .2)  170,  1*2 
170  KGO  s 3 

SUM  = SUM  - FI  • G1 
GO  TO  182 

172  SUM  = SUM  ♦ CMPLXCO.,  3.  / DELT6)  • (G1  • (FI  - .5  • (FI  • CMPLX(2 
1.,  3.  • OELTG)  - 4.  • F2  • CHPLXd,,  OELTG)  ♦ F3  • CMPLX(2.,  OELTG 
2))  / OELTG  * • 2)  - G2  • (F2  - ,5  • (FI  • CMPLX(2.,  OELTG)  - 4.  • 
3F2  ♦ F3  • CMPLX(2.,  - OELTG))  / OaTG  • ♦ 2)) 

182  OLOR  = R 
115  CONTINUE 

GO  TO  (180,  103) , IGO 

103  F3  * F2 
G3  = G2 
F2  = FI 
G2  = G1 

CALL  GEOHd,  I - 1,  2,  H(L)) 

G1  = CMPLX  (COS (NAVE  • R) , - SIN(MAVE  • R) ) 

Ft  X T(I  - 2)  • T(l)  * U • OEN  • M<I  - 1) 

OELTG  = HAVE  • (R  - OLOR) 

IF  (ABS(DELTG)  .GT.  .2)  184,  185 

104  SUM  * SUM  ♦ CMPLX(0.,  3.  f OELTG)  • (G1  • (FI  - .5  • (FI  • CMPLX(2 
1.,  3.  • OELTG)  - 4.  • F2  * CHPLXd.,  OELTG)  ♦ F3  • CMPLX(2,,  OELTG 
2))  / OELTG  • • 2)  - G2  • (F2  - .5  • (FI  • CMPIX(2.,  OELTG)  - 4.  • 
3F2  ♦ F3  • CMPLX(2.,  - OELTG))  / OELTG  • • 2)) 

GO  TO  188 

185  SUM  X SUM  ♦ 1.25  • FI  • G1  ♦ 2.  • F2  • G2  - .25  » F3  • G3 
188  00  190  K X 1.  4 

CALL  GEOMd,  1000  ^ K,  3,  H(L)) 

SUM  X SUM  ♦ 3.  • EH(K)  • ME(K)  • CMPLX (COS (MA VE  • R) , - SIN(HAVE  • 
1 R))  • OEN  • U / (SQRT(XS(K))  ♦ T(l)) 

190  CONTINUE 

SUN  * .33333333  • OELTX  • SUM 
0 X TX  / T(I  - 1) 

MM  X d.  - (^LX(Q,  Q)  • SUM)  • .5 
193  CONTINUE 

CALL  INOF  (H(L),X(I) ,NUT,FINOH) 

TS=NM*FINOH 

EdS,  L)  X CABS(TS)  * 1.257  * FREQ  t RO 
FdS,  L)  X - (CANG(TS)  - NAVE  • (RO  - X(I))) 

195  CONTINUE 
110  CONTINUE 

OISS  (IS)  X .001  * Xd) 

200  CONTINUE 

105  CONTINUE 


IF 

(LTOP  .LT. 

1)  GO 

TO 

120 

00 

250  L X 1, 

LTOP 

IF 

(ITCP  .LT. 

1)  GO 

TO 

125 

00 

240  IS  X 1, 

ITOP 
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FM  = FdS,  L)  • AUrCRO 
IFIIS.N6.NPALT)  GO  TO  9877 
AlTTNSVJITCT.L) =FM 
9877  CONTINUF 
?40  CONTINUE 
125  CONTINUE 
250  CONTINUE 
120  CONTINUE 
280  CONTINUE 
RETURN 
END 


SUBROUTINE  GROUNOCI,  K,  IGO,  HH) 

COMMON  /GRGUNO/ZDUM,  It  ZPDUM,  7P,  XDUH,  Xf  RO , DEN,  R,  U,  OMX,  MA 
Hi,  HZ,  2PP,  DELTA,  OELTAl,  OELTA2 
COMMON  / ZOTA/ARRAY 

COMMON  /S0i:0I/S(999)  , ORC999),  01 1999)  , LLH 

DIMENSION  ZOUM(6),  Z(1009),  ZP0UM(6),  ZP(1009),  XOUM(6),  X(1009), 
lARRAV  (NO ) 

COMPLEX  DELTA,  DEN,  OELTAl,  0ELTA2 
EQUIVALENCE  (ARRAY(l),  FKHZ),  (ARRAY(13),  A) 

CALL  INT  (K,I) 

FREQ  = FKHZ  • l.E-3 

RETURN 

ENTRY  GEOM 

HIT  = Z(  I)  ♦ HH 

IF  (I  .EQ.  K)  GO  TO  2C 

r = (X(I ) - X(K) ) / A 

6S  X A ♦ Z(K) 

GX  = A ♦ HIT 

cr  = coscT) 

ST  = SINtT) 

01  = OMCOS(T» 

R2  X SORT(2.  • GS  • GX  • OT  ♦ (HIT  - ZCK))  • • 2) 

ftl  = SORT(2.  • A * GS  • OMCOSIX(K)  / A)  ♦ Z(K)  • • 2) 

Rx  PI  4 R2*  RO 

U = X(K)  * RO  • SQRTd  ♦ ZP(K)  * • 2)  / (R1  * R2  • X(I)) 

IF  (IGC  .LT.  3)  U = U • (X(I)  - X(K)) 

PO  = (A  • OT  ♦ Z(K)  - HIT  * CT  ♦ GX  • ZP(K)  • A • ST  / GS)  / R2 

XK  = X(K) 

IF  (XK  .GE.  S(LLM))  GO  TO  12 
LLL  = 1 

IF  (XK  .LT.  S( 1) ) GO  TO  10 
LLMO  X LLM  - 1 
IF  (LLMO  .LT.  1)  GO  TO  100 
00  13  LLL  X 1,  LLMO 

IF  (XK  .LT.  S(LLL  1)  .AND.  XK  .GE.  S(LLL))  GO  TO  10 
13  CONTINUE 
100  CONTINUE 
GO  TO  21 

10  DELTA  X CHPLX(DR(LLL)  , DI(LLL))  • SORT  (FREQ  / 0.1) 

GO  TO  21 

12  DELTA  X CMPLX(OR(LLM) , DI(LLM))  • SORT (FREQ  / 0.1) 

21  CONTINUE 

DELTA2  X CHPLXd.,  - 1.  / (HAVE  • R2))  • PD 
OELTAl  X delta  - 0ELTA2 
19  DEN  * OELTA2  ♦ DELTA 
RETURN 
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20  U = 0. 

R2  = GX  - « - Z(I) 

RETURN 

EMO 


SUBROUTINE  6ETELV<LATX,L0NX,AMP,FAZ,EL VTN) 
COHHON/GET£LV/L(120I  «IHP(120) 

COMMON/UNPACK/NPACX(60I$COHMON/FLAGS/IXPT,NFLAG 
OIHENSION  NSINOXfSSOS) 

DATA  KPSl/1/ 

INOOIH=5800 

IF(KPS1.  EQ.DCALL  OPENNS  ( 2,  HS  INOX  .INOOIH,  0)  SKPS1=0 
H11=3777B$N8=377B 

IF(LATX.LT.A80000  .OR.  LATX. GE. 54010 0)  GO  TO  190 
IFCLONX.LT. 080000  .OR.  LONX. GE. 14000 0)  GO  TO  190 
NFLAG=0$LATSEC=NOO(L ATX,10Q)$LONSEC=NOO(LONX,100) 
LATDM=LATX-LATSEC«LONON=LONX-LOIISECtlF  (L ATSEC.GE.3 0» GOTO40 
LATS=LATOMtLATN=LATOH*30$GO  TO  60 
40  LATS=LATON*30SLATN=LATOH4'100SLATNsHOO(LATOMflOOOa) 

IF  (LATH.  EO. 5900)  L ATN=LAT ON4^4100 
60  lF(LONSEC.GE.30)GOTO80SLONN=LONOHSLONE=LONON«30$GOTO100 
80  LONH=LONOH'»30SLONE  = LONDH'«’100 
LONH=HOO  (LONDH. 10000) 

IF(LQNH.  EO.5900)  LONE^LONOHi-AlO 0 
100  CONTINUE 

INDNE  = INDEX (LAIN,  LONE) 

IF(NFLAG.EQ.l)  GO  TO  120 
CALL  REAOHS(2,NPACKy60f INONE) 

CALL  UNPACK 
rME=L (IXPT) 

IF (NOD (LATN, 10000) .EQ.fl.OR.HOO (LO*C  , 20 000 ) .EQ.O) 60T01 20 
I5E  = L(IXPT-1) 

call  REAOHS(2«NPACK,60tINONE-l)tCALL  UNPACK 
INH=L(TXPT)fISH=L(IXPT-l)tGO  TO  180 
120  CONTINUE 

INDNN=INDEX(LATN,LONM) 

IF(NFLAG.EQ.l)  GO  TO  140 
CALL  READHS(2,NPACK,60|INONM) 

CALL  UNPACK 
INM=L (IXPT) 

140  CONTINUE 

INOSE=INDEX(L ATS (LONE) 

IF(NFLAG.EO.l)  GO  TO  160 
CALL  REAONS(2|NPACK,60(INDSE) 

CALL  UNPACK 
ISE=L(IXPT) 

160  CONTINUE 

INOSMxINCEX(LATS.LONH) 

IF(NFLAG.EC.l)  GO  TO  180 
CALL  REAONS(2f NPACK,60fINOSH) 

CALL  UNPACK 
ISM=L  (IXPT) 

180  CONTINUE 

AHPxFLCAT(SHIFT(IHP(IXPT) ,-11) .A.H8)«. 001 
FAZ=FLOAT  (IHP(IXPT) .A. Nil) *.0  01 
IF(NFLAG.EO.l)  PRINT  3100 
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ELEAST«FLOAT(ISE)4R.OAT{(lNE-tSE)«(LATSEC-(LATSEC/30)«30>)/30. 

ELMEST>n.O«T(ZSM)«rLOATI(INM>ISM)«a«TSEC-(LATSEC/30)«3Q))/30. 

ELWTN>ELMEST»(ELEAST'>ELMEST>«FLOATILOMSEC-(LONSEC/30)«30)/3fl. 

RETURN 

191  PRINT  2970 
RETURN 

2971  FORHATf*  •,<iH****,*LAT,tON  REQUESTED  ARE  OUTSIDE  NAP  REGION  - NO  E 
ILEVATION  RETURNED**//) 

3199  FORHATC*  •,AH****,* INDEX  REQUESTED  EXCEEDS  SCANLINES  GENERATED,  TH 
lUS  PREVIOUS  ELEVATIONS  NILL  BE  USED*) 

END 


FUNCTION  INDEX fLAT*L ON) 

COHNON/FLAGS/IXPTfNFLAG 
LATS«NOO  ILAT,10Q) 

LATH*NOO(LAT,10090)-LATS 

LATOs(LAT-LATN-LATS)/10000 

LATN-LATN/100 

LONS>NOO(LON*100) 

LONNsNOO  (LON* 10 090) -LONS 
LOND> (LO  A-LOHN-LONS) /lOOO  9 
LONN>LONH/100 
ISUBLN«(LONO*5) 

ISUBLN>(NOO(ISUBLN*10)/S)*120 

ILONO>((LQNO-6)/2)*l<»A9*ISUBLN«LONH*2*LONS/3B*l 

INDEX>ILONO«(LATO-Aa)*299 

IXPT*LATN*2H.ATS/30*1 

IF  Cl NOEX. GT.5760 ) NFL  AG*1 

RETURN 

END 


SUBROUTINE  UNPACK 
CONNON/6ETELV/L ( 120 ) , INP ( 120 ) 

C ONNON/UNPACK/NPACK (6  0) 

DATA  NSKl/3777B/,NSK2/ir77777B/ 

00  20  NA*1*60$N-2*(NA-1)>1 
L ( N)  > AND  (NSKl, SHIFT  CNPACK  (NA)  *-4i9) ) *9 
INP(N)sANO(WR2*SHIFT(NPACK(NA),30)  ) 

LCN^DsANO (HSK1*SHIFT(NPACK1NA)  ,-19))*8 
IF(L(N).6T.290  0)L(N)>96SIFCL(N'»1).GT«2090)LCN*1)=96 
INPiH«l)sANO(NSK2,NPACK(NA)) 

20  CONTINUESENO 
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SUBROUTINE  GEOOKBl*  ALl,  AZIZ.  S.  BZ«  ALZ,  AZZl) 

CALCULATES  INVERSE  COHPUTATION  FORH  CSOOAMO,  1965) 

Bl.ALl.BZ.ALZ  ARE  GIVEN  DATA 

AZIZ.S.  AZZl  ARE  CALCULATED  DATA 

ai  = GEODETIC  LATITUDE.  RADIANS,  OF  FIRST  POINT 

All  > GEODETIC  LONGITUDE,  RADIANS,  OF  FIRST  POINT 

BZ  = GEODETIC  LATITUDE,  RADIANS,  OF  SECOND  POINT 

ALZ  « GEODEnC  LONGITUDE,  RADIANS,  OF  SECOND  POINT 

——AZIZ  * FORNARO  AZINUTH,  RADIANS,  FRON  FIRST  POINT 
—AZZl  > BACK  AZINUTH, RADIANS,  FRON  SECOND  POINT 

S = DISTANCE  (LENGTH  OF  GEODETIC  LINE)  BETHEEN  POINTS,  NETERS 

—SOUTHERN  LAHTUOES  AND  WESTERN  LONGITUDES  ARE  NEGATIVE 

COHHON  ZGENERAL/TPI,  PI,  PIZ,  PIN,  AO,  BO,  FL,  ESQ,  IFLAG 
DATA  (IFLAG  c 1) 

IF  (IFLAG  .EQ.  1)  CALL  SETUP 
IF  (A6S(B1I  .GT.  PIN)  GO  TO  1 
TBETl  = TAN(Bl)  • (1.  - FL) 

BETl  < ATAN(TBET1) 

GO  TO  Z 

1 CBETl  * 1.0  Z TAN(81)  / (1.  - FL) 

BETl  s ATANd.  / CBETl) 

Z CONTINUE 

ALLl  = ALZ  - ALl 
IF  (ALZ  - ALl  .EQ.  0.)  6,  9 
0 ALLZ  - 0. 

GO  TO  3 
9 CONTINUE 

f ALLZ  = ALZ  - ALl  - SIGN  (TPI,  ALZ  - ALl) 

• 3 IF  (A6S(ALL1)  .GT.  ABS(ALLZ))  5,  6 

I 5 ALL  > ALLZ 

t GO  TO  7 

I 6 ALL  > ALLl 

; 7 CONTINUE 

IZ  IF  (ABS(ALL)  .EQ.  0.  .OR.  ABS(ALL)  .EQ.  PI  .OR.  A8S(ALL)  .EQ.  TPI) 
I 1 10,  11 

10  ALL  * ABS(ALL) 

' 11  CONTINUE 

IF  (ABS(6Z)  .GT.  PIN)  GO  TO  16 
TBETZ  = TAN(BZ)  • (1,  - FL) 

BETZ  = ATAN(TBETZ) 

GO  TO  17 

16  CBETZ  s l.I  Z TAN(BZ)  Z (1.  - FL) 

BETZ  = ATANd.  / CBETZ) 

17  CONTINUE 

CBETl  = COS(BETl) 

SBETl  > SIN(BETl) 

CBETZ  = COS (BETZ) 

SBETZ  = SIN(BETZ) 

A = SBETl  • SBETZ 
B = CBETl  • CBETZ 
AB  = SBETl  • CBETZ 
BA  = SBETZ  * CBETl 
COSL  = COS (ALL) 

SINL  - SIN  (ALL) 

CPHI  = A ♦ B • COSL 

SPHI  * SORT  ((SINL  • CBETZ)*PZ  ♦ (BA  - AB  • COSL)**Z) 

IF  (SPHI  .GE.  C.  .AND.  CPHI  .GE.  0.)  ZO,  Z1 
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20  PHI  s ASINCSPHI) 

IF  tSPHI  .6T.  CPHI)  PHI  = ACOS<CPHI) 

GO  TO  30 

21  IF  (SPHI  .GE.  0.  .AND.  CPHI  .LE.  0.)  22«  23 

22  PHI  « PI  • ASIN(SPHI) 

IF  (SPHI  .CT.  ABS(CPHI))  PHI  = ACOS(CPHI) 

GO  TO  30 

23  IF  (SPHI  .LT.  0.  .AMO.  CPHI  .LT.  0.)  2A,  25 

24  PHI  < PI  - ASIN(SPHX) 

IF  (ABS(SPHn  .6T.  ABS(CPHI))  PHI  = TPI  - ACOS(CPHI) 

GO  TO  30 

25  IF  (SPHI  .LT.  0.  .AND.  CPHI  .6E.  0.)  26.  27 

26  PHI  « TPI  * ASIN(SPHI) 

IF  (ABSCSPHU  .6T.  CPHI)  PHI  = TPI  - ACOS(CPHI) 

60  TO  30 

27  CALL  EXIT 
30  CONTINUE 

C = a • SIHL  / SPHI 
FL2  = FL  • FL 
CONI  = FL  ♦ FL2 
COM2  = 0.5  • FU2 
C0N3  = SPHI  • CPHI 

C0N4  = PHI**2  / SPHI 

CONS  = C0N4  • CPHT 

EH  * 1.  - C • • 2 

RATIOl  = (1.  ♦ CONI)  • PHI  ♦ A • (CONI  • SPHI  - CON2  • COH4) 

1 ♦ EH  • (-0.5  * CONI  • (PHI  ♦ C0N3)  ♦ CON2  • C0H5)  - A • A • C0N2 

2 • C0N3  ♦ EH  • EH  • CON2  • (0.125  • (PHI  ♦ COH3)  - CON5  - 0.25  • 

3 CON3  • CPHI**2)  ♦ A • EM  • CON2  • (C0N4  ♦ CON3  • CPHI) 

S = RATIOl  • BO 

IF  (S  .LE.  l.E-41  S = 0. 

RATI  02  * com  • PHI  - A • CON2  • (SPHI  ♦ 2.  • CON4) 

1 ♦ 0.5  • EH  • COM2  • (-5.0  • PHI  ♦ CON3  ♦ 4.0  • CONS) 

ALAM  * RAT 102  • C ♦ ALL 
SALAH  3 SIN(ALAH) 

CALAH  = COS(ALAN) 

CTAZ12  s eA  - CALAH  * AB 
CTAZZl  = BA  • CALAH  - AB 
IF  (ALl  - AL2  .EO.  0.)  35.  39 

35  AZ12  * 0. 

AZ21  0. 

GO  TO  34 

39  CTAZ12  = CTAZ12  t (SALAH  • CBET2) 

IF  (CTAZ12  .EQ.  0.)  54.  55 

54  AZ12  = PI2 
GO  TO  56 

55  CONTINUE 

AZ12  = ATANd.  / CTAZ12) 

56  CONTINUE 

CTAZ21  = CTAZ21  / (SALAH  • CBETl) 

IF  (CTAZ21  .EQ.  0.)  57.  SO 

57  AZ21  = PI2 
GO  TO  34 

58  CONTINUE 

AZ21  > ATANd.  / CTAZ21) 

34  CONTINUE 

IF  (ALL  .GE.  0.  .AND.  CTAZ12  .GE.  0.)  40.  41 


39 


o o 


CTAZIZ  .LT 


1.)  4Z,  4»3 


i»0  AZIZ  < AZIZ 
GO  TO  50 

41  IF  (ALL  .6E.  0.  .AND. 
4Z  AZIZ  - PI  ♦ AZIZ 
60  TO  50 
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IF  (ALL  .LT.  1.  .AND. 

CTAZIZ 

• 

UJ 

• 

0.) 

44, 

45 

44 

AZIZ  X PI  . AZIZ 

GO  TO  50 

45 

AZIZ  X TPI  ♦ AZIZ 

50 

IF  (ALL  .GE.  0.  .AND. 

CTAZZl 

• 

111 

• 

0.) 

46, 

47 

46 

AZZl  X PI  . AZZl 

GO  TO  51 

47 

IF  (ALL  .GE.  0.  .AND. 

CTAZ21 

.LT. 

0.) 

40, 

49 

48 

AZZl  X TPI  ♦ AZZl 

GO  TO  51 

49 

IF  (ALL  .LT.  0.  .AND. 

CTA  Z21 

0.} 

sz. 

53 

5Z  AZZl  X AZZl 
GO  TO  51 

53  AZZl  * PI  ♦ AZZl 
51  CONTINUE 

AZIZ  > AHOOIAZlZt  TPI) 

IF  (AZIZ  .LT.  0.)  AZIZ  > AZIZ  ♦ TPI 
AZZl  > ANOO(AZZly  TPI) 

IF  (AZZl  .LT.  0.)  AZZl  > AZZl  *■  TPI 

RETURN 

END 


SUBROUTINE  SETUP 

COHHON  ZCENERAL/TPI,  PI,  PIZ,  PZ4,  AO,  BO,  FL,  ESQ,  IFLA6 
DATA  (TPI>6.Z831S530717959} 

PI  « TPI  / Z.0 
PIZ  » TPI  / 4.0 
PI4  X TPI  / 0.0 

C -AOxSENIMAJOR  AXIS  OF  CLARKE  SPHEROTD  OF  1066,  METERS 

C BOxSEHIHINOR  AXIS  OF  CLARKE  SPHEROID  OF  1066,  METERS 

AO  X 6.37eZ064E46 
BO  X 6.356Se30E«6 
— — FLxSPHEROIOAL  FLATTENING 

ESQxSECOND  ECCENTRICITY  SQUARED 

FL  x 1.0  - 00  / AO 

ESQ  X (AO  • * Z - BO  * * 2)  / BO  ♦ • 2 

IFLAG  X 0 

RETURN 

END 
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FUNCTION  OHCOS(X) 

C ONOCOS(X)=  I - COS<X) 

IF  (ABSIX)  .GT.  .15)  GO  TO  40 

IF  (X  .EO«  0.)  GO  TO  50 

S = X • X 

OHCOS  s T = ,5  • S 

R = 4. 

10  T=-T*S/(R*CR-  l.)> 

ONCOS  = CMCOS  ♦ T 

IF  (T  / ONCOS  .LE.  .5E-9)  GO  TO  51 
R = R ♦ 2. 

GO  TO  10 

40  ONCOS  = 1.  - COS(X) 

RETURN 

50  ONCOS  = 0. 

51  RETURN 
END 


SUBROUTINE  INT  (I,K) 

C I = POSITION  IN  X ANO  Z ARRAYS  ON  WHICH  TO  CENTER  CAtCULATIONS 

C K = POSITION  IN  ARRAYS  Z AND  ZP  TO  STORE  CALCULATED  YALUES 

CONNON  /GROUND/  ZOUN , Z .ZPOUN .ZP* XOUH,X ,R0 >DEN,R.U, ONX , MA VE,R2« ZPP 
l.D, ODDI, 0002 
CDHPLEX  DEN.Dt 0001,0002 

DIMENSION  Z0UH(6) ,ZC1009) ,ZP0UH(6), ZP(1009) (XODM (6 ) • X ( 1009) 

IM0=I-1 

IP0=I+1 

C*({Z(IPO)-Z(IMO))-IZ(I)-Z<IMO))*<X(TPO)-xaMO))/(X(I)-KlIMO)))/ 

1 ( (X(IPO)*X(IPO)-X(IMO)*X (IMO))-(X<I)*X(I)-X(IMO)*X(IMO))* 

2 (X(IPO)-X(IHO))/(X(I)-X(INO))) 
B=(tZ(I)-Z<IMO))-C*{X(I)*X(I)-X(IMO)*X (INO)))/(X (I)-X(IHO)) 
A»Z(I)-X(I)*<B<-C*X<I)) 

Z(K)=iA*X  (K)*CB«-C*X(IC)) 

ZP(K)=B+2.*C*X(X) 

ZPP=2.*C 

RETURN 

END 


SUBROUTINE  CNEVKEN(A,  FA,  NA,  X,  FX,  NX,  NPT,  KASE ) 

COMPLEX  FA,  FX,  FUNCT,  POLY 

DIMENSION  A(NA),  FA(NA),  X(NX)  , FX(NX)  , FUNCTdS),  ABSCdS),  0IF(1 
15),  P0LY115) 

1 FORMAT  (1H1/50H  THE  X VALUES  ARE  NOT  ARRANGED  IN  ASCENDING  ORDER./ 
1/5X,4HI  = E20.9,5X,7MX(I)  » E2B.9,5X,4HJ  = E20.9,5X,7HA t J)  = E20.9 
2///14X,lHX,2DX,4HF(X)) 

2 FORMAT  (5X,2E20.9) 

200  FORMAT  (IHl*  THERE  ARE  NOT  ENOUGH  POINTS  IN  THE  GIVEN  ARRAY.*) 

205  FORMAT  (IHl*  THE  GIVEN  VALUE  OF  NPT  WAS  GREATER  THAN  NA,*/*  THE  SO 
IBROUTINE  REDUCED  NPT  TO  WITHIN  THE  LIMITS  OF  NA.*) 

210  FORMAT  (IHl*  NPT  WAS  INITIALLY  GREATER  THAN  15  - NPT  «I5,*/*  THE  S 
lUBROUTINE  SET  NPT  TO  15  ANO  CONTINUED.*) 

LOOP  = 1 

IF  (NPT  - 15)  3,  3,  9 

9 PRINT  210,  NPT 
NPI  « 15 
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3 ir  (NPl  - NM  8,  8t  *» 

% MPf  * NPT  - Z 

IF  CNPT  - 1)  5f  5,  6 

5 PRINT  200 
CALL  EXIT 

6 GO  TO  (7,  3)(  LOOP 

7 LOOP  » 2 
PRINT  20G 
GO  TO  3 

8 NPT2  * NPT  / 2 

GO  TO  (11,  12),  KASE 

11  NSTART  * 1 
NX  > 1 

GO  TO  16 

12  NSTART  - NA  - 3 

TEST  A (NSTART)  *■  (A(NSTART  ♦ 1)  - A(NSTART))  / 2 
IF  (NX  .LT.  1)  GO  TO  170 
DO  16  I > 1,  NX 
IF  (X(T)  - TEST)  16,  16,  13 

13  NX  > I 
GO  TO  16 

16  CONTINUE 
170  CONTINUE 
16  NSTOP  « NA  • 1 

IF  (NX  .LT.  NX)  GO  TO  17S 
DO  125  I 3 MX,  NX 
IF  (X(I)  - A(l))  135,  15,  10 

10  IF  (X(I)  - A(NA))  25,  20,  130 

15  FX(I)  * FA(1) 

GO  TO  125 
20  FX(I)  = FA(NA) 

SO  TO  125 

25  IF  (NSTOP  .LT.  NSTART.)  GO  TO  18R 
DO  as  J > NSTART,  NSTOP 
IF  (X(I)  - A(J))  32,  35,  30 

30  IF  (X(I)  > A(J  *■  1))  65,  60,  85 

32  II  ■ I - 1 

PRINT  1,  1,  X(I),  J,  A(J) 

IF  (ID  36,  36,  33 

33  PRINT  2,  (X(N),  FX(N),  N > 1,  ID 
36  CALL  EXIT 

35  FX(I)  « FA(J) 

NSTART  s J 
GO  TO  125 

60  FX(I)  a FA(J  ♦ 1) 

NSTART  s J ♦ 1 
GO  TO  125 
65  NSTART  > J 

IF  (ABS(X(I)  > A(J))  - ASS(X(I)  - A(J  * 1)))  50,  50,  55 
50  JJ  « J 
GO  TO  60 
55  JJ  ■ J ♦ 1 
GO  TO  60 
85  CONTINUE 
180  CONTINUE 

60  IF  (JJ  - NPT2)  135,  135,  70 
70  IF  (JJ  * NPT2  * NA)  00,  80,  130 
80  KX  a JJ  - NPT2  - 1 
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90  IF  (NPT  .LT.  1)  GO  TO  185 
DO  95  K s 1,  NPT 
KK  * KK  ♦ 1 
FUNCT(K)  » FA(KK) 

ABSC(K)  s A(KK) 

95  OIF(K)  « AeSC(K)  • X(I) 

185  CONTINUE 

NTOP  » NPT  - 1 
LL  * 1 

100  IF  (NTOP  .LT.  1>  GO  TO  190 
DO  105  L = 1,  NTOP 
LLL  > L ♦ LL 

105  POLTILl  = (FUNCT(L)  • OIF(LLL)  - FUNCT  (L  ♦ 1>  • OIF(Ln  / 
ID  - ABSC(D) 

190  CONTINUE 

IF  (NTOP  - 1)  120.  120t  110 
110  00  115  N * 1,  NTOP 

115  FUNCT (HI  > POLY(N> 

NTOP  a NTOP  - 1 
LL  » LL  ♦ 1 
GO  TO  100 
130  INC  * - 1 

KK  = NA  « 1 
GO  TO  140 
135  INC  z 1 
KK  > 0 

140  IF  (NPT  .LT.  1)  GO  TO  215 
00  145  K s If  NPT 
KK  s KK  « INC 
FUNCKK)  s FA(KK) 

ABSC(K)  s A(KK) 

145  OIF(K)  X AeSC(K)  - X(I> 

215  CONTINUE 

NTOP  X NPT  - 1 
LL  = 1 

150  IF  (NTOP  .LT.  II  60  TO  220 
00  155  L X 1,  NTOP 
LLL  X L ♦ LL 

155  POLY(L)  X (FUNCT(I)  • OIF (LLL)  - FUNCT (L  ♦ 1)  • OIF(LL)) 
ILL)  - ABSC(LL)) 

220  CONTINUE 

IF  (NTCP  - 1)  120f  120.  160 
160  00  165  N X 1,  NTOP 

165  FUNCT(M)  * POLY(M) 

NTOP  X NTOP  - 1 
LL  X LL  ♦ 1 
GO  TO  150 

120  FX(I)  X POLY(l) 

125  CONTINUE 
175  CONTINUE 
RETURN 
END 


(ABSC(LL 


/ (ABSC(L 
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COMPLEX  FUMCTION  FLEAFCMAVE,  HI,  H2,  XD,  OELTAR) 

COMPLEX  TEMP,  Q,  I,  Z2,  ZZ,  HMERF,  MERFZ,  HERF,  ZMERF,  OELTAR 
HO  = H2  - HI 

TEMP  * (0.7iZliJ6Z812,  - 0.7071067812)  • SQRT(.5  • MAWE) 

X02  = SQRT(XO) 

Q = - TEMP  • HO  / X02 

Z = TEMP  • OELTAR  * X02  ♦ Q 

ZZ  * - Z 

ZI  > AIMAG(ZZ) 

IF  (Zl  .LT.  0.  .OR.  (A8S(REAL(ZT) > .LT.  6.  .ANO.  ZI  .LT.  6.1)  GO  T 
10  10 

Z2  = ZZ  • • 2 

HMERF  = (Z2  - 2.1  / (ZZ  • (Z2  - 3.5)) 

GO  TO  12 

C HERF  - COMPLEMENTARY  ERROR  FUNCTION 

11  HERFZ  = HERF(ZZ) 

HHERF  = ZZ  - 0.5  • MERFZ  / (ZZ  • HERFZ  ♦ (0.,  - 0.56A18950)) 

12  ZHERF  * Z ♦ HHERF 

FLEAF  = (Q  • ZHERF  - 0.5)  / (Z  • ZHERF  - 0.5) 

RETURN 

END 


COMPLEX  FUNCTION  HERF(ZZZ) 

COMPLEX  Z,  ZZZ,  ZV,  V,  Z2 , C,  H,  S 
DIMENSION  C(12),  H(5,  4) 

EQUIVALENCE  (S,  C(12)) 

DATA  (C(l)  = (.0,  - .5641095835)) 

DATA  H/(1.,0.0), 

1 (3.678794411714423E-01,6,071577058413937E-01) , 

2 (1.831563ee8073410E-02,3,400262170660662E-01) , 

3 ( 1.2340 9e04BS66788E-04, 2. 011573170 376Q04E-I1 ) , 

4 (1.1253Sl747192646E-07,1.45953589900152eE-01) , 

5 (4.275835761558070E-01,0«000O0OO000Q0000E«00)  , 

6 (3.047442052569126E-01,2.oa21893e2a20316E-01) • 

7 (1.4023«58136b2779E-01,2.222134401798991E-01) , 

8 (6.531777728904697E-02,1.739183154163490E-01) , 

9 ( 3. 6 2814Se46998864E>a 2,1. 358389510 006551E-01 ) , 
A ( 2. 55395676310 5058E-a 1,0. 0000 OOOCO 0 00  0 00 E« 00 ) , 
B (2.1849  <6152748907E-01, 9. 2997809392601862-02) , 
C (1.47952759512015BE-01,1.311797170842178E-01) , 

0 (9.2710766426443322-02, 1. 283169622 282615E-01) , 
E (5. 9686929610445902-02, 1.1321005612448822-01), 
F (1.790011511813930E-01,O.OOOOOQOOOOOOOOOEtOO) , 
G (l,642611363929861E-01,5.019713513524966E-02) , 
H (1. 307574696698522E-01, 8. 1112650477454722-02) , 

1 (9,6402505503044392-02,9.1236326004212582-02) , 
J (6. 979096164964750E-02,8.934000024036461E-02)/ 

XX  = REAL(ZZZ) 

YY  = AIHAGdZZ) 

X = ABS(XX) 

Y = AeS(YY) 

Z = CMPLX(X,  Y) 

LZ2  = 0 

IF  (X  .GE.  4.5  .OR.  Y .GE.  3.5)  GO  TO  100 
I = X ♦ .5 
J = Y ♦ .5 

V = CMPLX(FLOAT(I) , FLOAT(J)) 
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ZV  * Z - V 

C(2)  = H(I  ♦ 1,  J ♦ 1) 

AT  s 0. 

03  10  I * 3>  12 
AI  = AI  - .5 

C(I)  * (W  • C(I  - 1)  ♦ Cd  - 2))  / AI 

10  CONTINUE 
J = 12 

03  11  I = 2,  11 
J = J - 1 

11  S = S • ZV  ♦ C(J) 

20  IF  (YY  .GE.  0.)  GO  TO  30 
IF  (LZ2  .EQ.  01  Z2  = Z • Z 
S = 2.  • CEXP(  - Z2»  - S 
IF  CXX  .GT.  0.)  S = CONJGfS) 

GO  TO  200 

30  IF  (XX  .LT.  0.)  S = CONJG(S) 

200  HERF  3 s 
RETURN 
100  LZ2  = 1 

Z2  = Z • Z 

S = Z • ((0.,  0.4613135279)  / (Z2  - 0.1901635092)  * (0.,  0.0999921 
16168)  / (Z2  - 1.7844927485)  ♦ (O.t  0.0028630938748)  / (Z2  • 5.5253 
24374379) ) 

GO  TO  20 
END 


function  CANG(Z) 

COMPLEX  Z 

DATA  (PI=3. 14159265358979) ,(PIHA=1. 57079632679489) 

XsREALtZ) 

Y=AIMAG<Z) 

IF  (X)  2C>30tlO 
10  CANG=ATAN2(Y,X) 

RETURN 

20  IF  (Y.NE.O.)  GO  TO  10 
CANG=PI 
RETURN 

30  IF  (V.GT.O.)  GO  TO  40 

IF  (Y.LT.O.)  GO  TO  50 

CANG^O. 

RETURN  I 

4C  CANG-PIHA 

RETURN  j 

50  CANG=-PIHA  1 

RETURN 

END  { 

1 

I i 

I I 


45 


! 


L__ J 


I 


SUBROUTINE  INDF  (HriyXX,NUT,F) 

C THIS  SUBROUTINE  CALCULATES  THE  INDUCTION  FIELDS  FOR  E SUB  R ANO 
C H SUB  PHI 

C THESE  induction  FIELDS  ARE  FOR  POSITIVE  TINE  FUNCTION 

C NUT  > 1 GIVES  INDUCTION  FIELD  FOR  E SUB  K (TAU  BAR) 

C NUT  « 0 GIVES  INDUCTION  FIELD  FOR  H SUB  PHI  (LOOP) 

C FZ  s INDUCTION  FIELD  FOR  E SUB  R 

C FH  X INDUCTION  FIELD  FOR  H SUB  PHI 

CDNHON/INDUCT/HAVE 
DOUBLE  THETAO,SINTHtCOSrH,R,CONS 
CQHPLEX  GZ«FZfFH>F 
36739  £6 
TPIs6«2eziaS307 
C-2.997925ES 
IF(XX.LE«S.)3,4 
3 PRINT  5 

5 FORHAT  (//•  IN  INDF.  DISTANCE  IS  ZERO  OR  NEGATIVEt  XX  > •,E20. 10) 
CALL  EXIT 
V THETAD>1.DB«XX/A 
SINTHxOSIN(THETAO) 

CaSTHsOCOS(THETAD) 

R^A^HH 

C0NS*A*R*SINTH**2 
D2»R*R^A*A-2.*A*R*COSTH 
IF(02.GT.O.)  GO  TO  30 
ODsXX 
02*XX*XX 
GO  TO  AO 
30  CONTINUE 
00=SaRT(02) 

Al  CONTINUE 
03s00*02 
DVs02*02 

FZR=-2.*C0STH/0O»3.*CONS/O3 

FZI*(NAVE*CONS*2.*COSTH/NAVE)/02-3.*CONS/ (OA'NAVE) 
FZ=CHPLX(FZR,FZI) 

GZcCHPLX(0.,NAVE) 

FZ=FZ/GZ 

FHR*R*SI NT ►•00703 

FHI*NAVE*R*SINTH*00/02 

FH*CNPLX(FHI.>FHR) 

C1=2.E-7*TPI*NAVE*C 

FH=FH/C1 

IF(NUT)t,2 

1 F«FZ 
RETURN 

2 F*FH 
RETURN 
END 


3.7  Sample  Test  Run 

A typical  data  output  listing  is  illustrated  below.  The  transmitter  coordinates 
which  are  read  in  from  the  data  card  deck  are  printed  for  verification  purposes. 

The  geographic  coordinates  of  the  desired  target  or  receiver  are  next  listed.  Array 
data  used  in  the  calculations  are  then  printed.  TOA  is  the  travel  time  from  trans- 
mitter to  receiver  in  psec  and  recorded  in  order  of  master,  slave  1 and  slave  2 to 
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target  respectively.  TDl  and  TD2  are  the  final  desired  LORAN  coordinates  in 
tisec.  EDI  and  ED2  are  the  slave  emission  delays  in  jisec  obtained  from  the  trans- 
mitter card  input  data.  Additional  output  data  recorded  as  TOA  above  include: 

DISTDUM  - Calculation  of  distance  by  Hufford's  technique  in  KM. 

TIMDUM  - Calculation  of  secondary  phase  factor  in  psec. 

DISTSOD  - Calculation  of  distance  by  Sodano  technique  in  KM. 

TPW  - Travel  time  of  primary  wave  in  jisec. 


3.8  Output  Listing 


T <LE  '-nnenTMAT^fc; 

LCLCJTIX-M UlXIXUilE LO'jai  T UDi. 

(nEG-MTM-crr., 


Mfl'tTVO 

-SLU'tEJ. 


prrcT|/cp  cCO^nTNarrc 

...  -LOCAITCK L-fiXIIUO£ LOflGl.TUO£ 

tPrc-uitg-crr)  OTr,-1TN-'E') 

REC'^IVFP  E ‘^r  31  29.1’6N  » -'7  u.irjE 


.ICCa"  fl2Et33 
.lC'’D''?<tE  + 01 
-«-3.3El-2  78E.*llZ 

2i.c  + nn 
♦ ';ic'L-a2E*3i 

.nno"  I'T*-''! 

. n 3c«.n3 

. m n TC  t-Pli 

.icco'm'^+ra 
.TtinanE-iQit 
. 1 cr"'  3"E*'’i, 
.ascciarEtcs 
.7liqi'3l,7E  + f'7 
T 1 nfinc  ff  nF*ni 


tHAHl 


TOA 

- r7i;  7C70=:fiAAZa.^  + "T^  . ^ 1 «-»A  r TAe;r  4.r  ^ ^ 

- i 1 •*  , 

ini  . 

= . 1' 7Ur  7777c;<tqf,»F+ '•S. 

Tnp 

= B1f,C7(171  7Firr*iq. 

EDI 

= *17  K37i.7F*7q^ 

tna 

-s-  .2'.1& itll-EtCE-. 

ARSaVllEi  FJ<a7..at-FP«^Qll^NClt  

4R0&y(21  = -la  = AfF  TNPFX  oc  a-FRaCTIOK  ST  GFCUND 

&fi-RA.y<31  = -aMfl.X  HaXTMilM  niSJAMC.f—'lEJltEFM.  IP  AMS.'tl  r T£9 

SRP.flVdi)  = PIMC  = '"ISTfiMCE  INCPEH^NT  IN  KM 
fl5g.aYiEA--=-  NSTaRT—.-IJiin’^X-  0F_  TUF-  FIRST  OISTaNPE  ai 

apRaytf.)  = mzinc  = '^he  field  is  to  te  fquno  as  a 

DiSJAtij::£.-.Aa£a.!(_--^-  nziuc-  must-  ‘iewer  nE 

iDoarc’)  = INDEX  OF  THt  lasT  PisiaNOE  ar  which  the 

ARRAYIFE  ZMTN-S  .1  HF  MINT  MU.1.ALJXTUIIE- ABDutE  THE 

*RRav(q»  = 7TMC  = the  aLTITUOE  TNC°EMENT,  METEFS 
ARRJiYllCl.  I 7MaY  I TM=-  »axTM|iH  -4LIJ.TU0£,_  .MEIE.RE 

ePRavtii)  ri  FLar  = the  FacTop  used  in  the  Ftai  EaciH 

apoaxti?)  = ai  PHa  = ucpTTrai  i apqr  Fartop 

apot.Ytl'')  = effective  EaPTHS  RaDIUS 
anpaxciLi  = miit 


L9  76  1«."1''N  ’ 1 9 36. ’’6E 

6 7 iq  l-».F.F7N »,  .-3  w".  E 

<*»  IE  <»?,9'’DN  11  77  i,a.2f3F 
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\.  DKSCRIPTION  OF  PART  II 


1. 1 Overview 

The  electrical  properties  of  the  earth  are  converted  into  a ground  surface 
impedance  in  PART  11  of  the  program  as  shown  in  Figure  7.  The  data  on  the  input 
tapes  from  PART  I are  loaded  on  to  a disc  so  that  at  each  index  point  information 
on  terrain,  geology,  and  soil  is  available.  This  information  is  employed  to  develop 
a 3 -layer  model  of  the  ground  as  shown  in  Figure  8.  The  time  of  arrival  of  the 
LORAN  pulse  at  the  receiver  is  controlled  by  the  electrical  characteristics  of  this 
model  through  a parameter  called  the  surface  impedance.  This  parameter  is  the 
ratio  of  the  horizontal  components  of  the  electrical  to  the  magnetic  field  components 

assuming  the  groimd  is  isotropic;  the  calculation  of  this  parameter  is  described  by 

T V,,  12,  13 

Johler.  ’ 

The  development  of  the  3 -layer  ground  model  requires  additional  investigation, 
in  particular,  the  constants  to  be  used  in  the  postulated  second  or  saturated  layer. 
The  electrical  characteristics  of  groimd  or  other  medium  may  be  expressed  by 
three  constants,  the  relative  permeability,  the  dielectric  constant,  and  the  conduc- 
tivity. The  relative  permeability  can  normally  be  regarded  as  unity.  From  the 

digitized  data  from  PART  I,  the  conductivity  of  the  top  unsaturated  ground  or  soil 

14  IS  16  IV  18 

and  basement  rock  can  be  determined.  » > > > Most  soil  maps  also  contain 

information  on  top  soil  depth  which  is  included  in  the  look  -up  table.  The  rock 
depth  is  assumed  to  be  infinite  since  the  100  kHz  wave  does  not  penetrate  below 
300  meters.  The  conductivity  and  depth  of  the  second  layer  are  somewhat  uncer- 
tain and  various  articles  in  the  journal  titled  Geoexploration^^'  have  been  used 
to  obtain  information  on  second  layer  parameters  as  a function  of  the  top  layer. 

The  values  listed  in  Figure  8 are  suitable  for  Europe.  In  Iran,  for  example, 
sigma  2 is  larger  than  sigma  1.  The  dielectric  constant  has  been  selected  as  15 

over  land  and  80  over  water.  An  empirical  relationship  between  conductivity  and 

21 

dielectric  constant  has  been  developed  by  Hanle,  and  may  also  be  used  in  the 


Due  to  the  large  number  of  references  on  this  page,  the  references  will  not  be 
footnoted.  See  References,  page  84. 


48 


r 


I 


Figure  7.  Data  Base  of  Ground  Electrical  Properties  - Part  II 


Zn 

ELEVATION  DATA  OFF  MAP 

dj  UNSATURATED  <^1  (MAP)  djtMAP) 

SOIL 

d2  SATURATED  <^2  ' ^2  “ ‘*2  ' ^1 


d3  '^3  (MAP)  ^3  ' 50  d3  ■ 00 

ROCK 

Figure  8.  Ground  Model  for  Surface  Impedance  Calculators 

program.  The  effects  of  temperature  changes,  vegetation,  buildings,  and  other 
objects  on  the  surface  of  the  earth  have  not  been  considered. 

The  input  data  to  PART  11  is  on  a standard  one-half  inch,  400-foot  magnetic 
tape  containing  information  on  elevation  in  meters,  bedrock,  and  soil  identification 
number.  The  input  format  for  the  data  is  as  follows: 

Elevation  - 11  bits  - ±7  meters 

Blank  - 3 bits 

Rock  Number  - 8 bits  - 256  types 

Blank  - 2 bits 

Soil  Number  - 6 bits  - 64  types 


49 


i 

j 

This  totals  to  a 30-bit  data  word.  Each  CDC  word  contains  60  bits,  therefore, 
two  data  words  or  points  can  be  accommodated  by  a single  CDC  word.  For  30-bit 
word  computers  (IBM  3 60,  UNIVAC,  Honeywell),  one  data  point  will  be  accommo- 
dated by  each  word  allowing  flexibility  between  machines.  There  are  60  words  in 
each  record  or  2 X 60  - 120  data  points.  With  a 30  arc  second  grid,  each  record 
covers  one  degree  in  latitude  of  the  grid. 

The  data  is  recorded  on  to  the  tape  in  the  following  sequence:  .Starting  at 
6°  E Long,  48°  N Lat,  in  the  southwest  corner  and  go  north  for  one  degree  in 
latitude  (one  record),  then  repeat  for  an  equal  line  incremented  in  longitude, 
parallel  and  30  arc  seconds  to  the  east  of  the  first  line.  This  is  repeated  until  the 
last  longitude  line  of  the  two  degree  standard  longitude  sectors  is  reached.  For 
this  particular  problem,  this  would  mean  after  2 X 60  X 2 = 240  records.  The 
latitude  is  then  incremented  by  one  degree,  returned  to  the  original  longitude,  and 
the  entire  process  repeated  for  another  240  records.  At  the  east  edge  of  the  final 
latitude  sector  (53°  to  54°),  the  pointer  returns  to  8°  N Lat  and  begins  another  1° 
latitude  scan  from  8°  E to  10°  E longitude.  For  the  six  degree  latitude  and  eight 
degree  longitude  area  coverage  in  Germany,  there  will  be  a total  of  960  X 6 = 57  60 
records  which  contain  5760  X 120  - 691,  200  data  points. 

The  output  data  tape  of  PART  il  which  is  the  input  to  PART  III  has  the  following 
format; 

Elevation  in  meters  - ±7  meters  - 11  bits 

Impedance  Amplitude  - ±2  ohms  - 8 bits 

Impedance  Phase  - ±.001  radians  -11  bits 

The  elevation  data  is  passed  directly  from  input  to  output  tape.  The  output  imped- 
ance is  returned  to  the  disc  for  storage  in  the  same  location  as  the  original  data 
and  then  recorded  on  magnetic  tape  for  shipment  to  the  PART  III  computation  sites. 

1.2  Sultrouline  Drsoriplion 

4.  2.  1 PROGRAM  CONIMP 

(a)  Controls  impedance  program. 

(b)  Call  subroutine  SETUP. 

(c)  Reads  and  writes  on  disc  soil,  geology  and  elevation  identification  data 
off  input  tapes. 

(d)  Call  subroutine  REPK. 

(e)  Writes  impedance  data  on  disc  and  output  tape. 

4.2.2  SUBROUTINE  SETUP 

(a)  Sets  relative  permeability  of  three  layers  to  unity. 

(b)  Reads  and  stores  ground  conductivity  tables  from  input  cards  into  look-up 

file. 
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L 4.2.3  SUBROUTINE  KEPK 

I 

] (a)  Determine  parameters  for  3 -layer  ground  surface  impedance, 

j (1)  Dielectric  constant  equals  fifteen  overland  and  eightly  over 

[ (2)  Determine  depth  of  first  and  second  layer.  .Set  third  layer  to  infinity. 

! (3)  Set  sigma  2 as  a function  of  sigma  1. 

' (b)  Call  subroutine  GRIM  for  impedance  calculation. 

4.2.4  SUBROUTINE  GRIM 

(a)  Calculates  surface  impedance  for  3 -layer  ground  model. 

(b)  Returns  complex  impedance  (AMP,  FAZ)  to  subroutine  RFPK  via  argument 

list. 

t.:!  Kesislivily  Table.s 

' For  subroutine  SETUP,  a look-up  table  is  required  to  convert  type  of  soil  and 

rock  into  resistivity  values.  The  bibliography  at  the  end  of  this  section  lists  avail- 
able sources  for  these  values.  Table  1 is  a list  of  average  values  for  the  top  or 
unsaturated  soil.  An  approximate  value  is  usually  assigned  in  the  map  legend. 

An  inspection  of  the  geological  map  of  Europe  reveals  characteristic  rock  types 
which  include  anorthosite,  granite  gneiss,  gabbro,  Precambrian  igneous,  and 
metamorphic  rock.  The  associated  resistivities  for  various  rock  types  are  defined 
in  Table  2. 


Table  1.  Ground  Resistivity  Values  for  Typical  Soils 


Terrain 

Conductivity 
(mhos  -meter) 

Resistivity 

(meter-ohms) 

.Sea  water 

5 

0.  2 

Fresh  water 

8 X 10'^ 

123 

Dry,  sandy,  flat  coastal  land 

2 X 10"^ 

/ 

3oo 

Marshy,  forested  flat  land 

8 X 10'^ 

123 

Rich  agricultural  land,  low  hills 

1 X 10'^ 

1,  00,0 

Pastoral  land,  medium  hills  and 

5 X 10'^ 

200 

forestation 

Rocky  land,  steep  hills 

2 X 10“^ 

500 

Mountainous 

1 X 10'^ 

1,  000 

Cities,  residential  areas 

2 X 10'^ 

500 

Cities,  industrial  areas 

1 X 10'‘* 

1,  oooc 
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t.4  Pro)?ttin  l.istin):  fur  Part  II 


program  CONIMP 


i'  program  C0NIMP(TAPE3,  YAPE2,  INPUT=128,OUTPUt=i28)' 

_ COMMON/ THREE /F ,P HI , SI GMA 1 , E 21 , AMU  1 , Z 1 , S I GM A 2 , E E2 , AMU2 , Z2 , S IG M A3 , 

~ ' .E23,AMII3 

COMMON/COOE/SLCD  (65, 3)  ,G0  ,GLCO  (1<.«) 

5 DIMENSION  NUPK (6 B ) , NO XEQ ( 58 0 0 ) 

OAT  A MSi<6/Z7B/,MSKF  15/77/77000000000  0000008/ 

LNLM=g60};LATLM  = 6 
LNLM=2CILATLM=1 
NbXSIZ=LNLM*LATLMr-40 
_ _CALL  SETUP 

CALL  OPENMS(2,NOXEO  ,5800 ,0) 

LNLM=576n 

DO  100  ILN=1,LNLM 
NOXL=ILN 

15~  BUFFERIN(3,1)  (NUPK( 1)  ,NUPK(60I )tIF(UNIT  131) 10,20 , 30 

3^0  PRINTl,  ILN 

1 FORMAT! *0  PARITY  ERROR*, 19) 

10  CONTINUE 

IF(MOO(  ILN,  1440)  . LT.IO)  PRINT  9,  ILN,  NDX L , NUPK'  (1)  , NUPK  (2) 

20  _CALLJ^P<JNUPK,MSK^,MS)^F15)  

IF(MOO(  ILN,1440)'.LT.iO)  PRI N T 99  , N UP  K < 1)  ,NUPK(2) 

ICVER=1 

IOVER=C 

70  CALL  WPITMS(2,NUPK,60,NOXL ,IOVER) 

25  9 FORMAT(*0  •,2I8,2022> 

99  FORMAT! 2022) 

100  CONTINUE 
20  STOPJFND 


SUBROUTINE  REPK 


1 SUBROUTINE  REPK ( NUP K , MSK 6, mSKFI 5) 

C0MM0N/THREE/F,PHI,SIGMA1,E21,AMU1,Z1,SIGMA2,E22, AMU2, Z2,SIGMA3, 
.£23,AMii3 

C0MM0N/C0DE/SLC0!65,3) ,G0 ,GLCD(148)  

5 DIMENSION  NUPK(60), iABilZO) 

DATA  MSK8/377B/ 

MSK30=37770 0000 03777000 0003 

DO  ^ 1=1,60 

LLEV=NUPK (I)  . A,  MSK30JLEV=SHIFT !LLEV, 1) 

JJ LS2=NUPK( I) . A.MSK6?LG2=MSK8. A .SHIFT (NUPK!I) , -6) 

LS1=MSK6.  A.SHIFT!NUFK!I)  ,-3C)  J LG  1 = MS K8  . A . SHI  F T ( NUP K ! I ) , - 36 )" 

no  99  I2=l,2$IF!I2.EQ»2)LSl=LS2tIF(I2.EQ.2)LGl=LG2 

E23=E22=E21=15. 

IF!LS1.  EC.O .OR.LSl .EQ.22.0R .LSl.EO. 12.0R.LG1 .GE. 146) E2  2=E21=80. 

15  SIGMAl=SLCn!LSl*l,l) 

71=SLCr  !LS1*1,2)*0.75 

IF! Z1 .IT. 2.)  Zl=2. 

IP!71.lT.SLCO!LSl»l,3))Zl  = SLCO!LSl»l,  *) 

IF(STGMAl.EO..01)  E22=E21  = 80. 
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20  Sir,MA2=0.  5*S  ir.MAlCZ?=<..*  Z1 

IP  {f2?.E0.9'3.)5;iGMa2=SI&Mai 
IL4T=2» T+I2-2 

5IGMA3=GLCri(  LGl ) JCALL  GRI M ( AMP,  FAZ  ) I lA  MP=A  M P$  I FA  Z=F  A Z 
N1=30* (2-12) 

2G  99  L EV=OR(  SHIFT  (I  AMP,  NB-fll)  .SHIFT  ( IFAZ,  NB)  , LEVTiNUPKt  f)  =LE\/ 
100  CONTTNUPJENT 


SUBRBUTIMF  GRIM 


1“  ■subroutine'  grim (AMPIMF,  FA^TMP)  “■ 

C F IS  rpFOJENCY  IN  HERTZ 

C OHI  IS  THE  angle  of  INCIORNCE  IN  RAOIANS 

C SIGMAl  IS  THE  CONDUCTIVITY  OF  THE  UNSATURATEO  SOIL. 

t;  C E?1  IS  THE  OIELECTRIC  CONSTANT  OF  UNSATURATEO  SOIL. 

C 71  IS  the  depth  in  METERS  OF  THE  UNSATURATEO  SOTL) 

C SIGMA?  IS  the  conductivity  OF  "THE  SAfURATEO  SOIL. 

C E22  is  the  OIELECTRIC  CONSTANT  OF  SATURATED  SOIL. 

C Z2  IS  the  DEPTH  IN  METERS  OF  THE  SATURATEO  SOIL. 

10  C SIGMAJ  IS  THE  CONDUCTIVITY  OF  PEOROCK. 

C E2T  IS  THE  OIELECTRIC  CONSTANT  OF  BEDROCK. 

C 2SOIL  IS  THE  TOTAL  DEPTH  OF  SOIL  ABOVE  BEORCCK  (Zl  ♦ Z2) 

COMPLEX  T2N,T2D 

rC“PLEX  AlO,  A20,  All,  A12,  A13,  A21,  A22,  A23,  A32,  A33, 

15  1A34,  AU2,  A1»T,  A44,  AKl,  AK2,  ETAl,  ETA2,  ETA12,  ETA2?, 

27ETA1,  7PTA12,  7ETA2,  7FTA22,  TEMP,  RENUM,  REOEM^RE,  PMI.BETAl 
COMPLEX  SAV1,SAV2,SAV3,SAV4,SAV5,SA V6 

COMPLEX  AK3,  FTAT,  FTAT2,  ZETA3,  7ER01,  A54,  A55,  A64, 

1A65,A66,  T2,  T3N,  T70,T3  , A35,  A45,  A55 

20  COMMON  /THREE/  F , PHI  j_SIGmA1,  E21,  AMU^,  71, _ ''IGMAZ, 


1E22,  AMU2,  7’,  SIGMA3,  E23,  AMU3 
data  (C  = 2 .B97B2SF8),  ( T P I = F . 2 83 1 8 53 0 7 1 796) 

OATA(  EO  = 6.‘85418F3367321E-12) 

DATAfPI  = 3.141592654)  ,{PI2  = -1.570796327) 

25  DATA  ISKIP/1/ 

IF  ( ISKIP. ED. 0) GO  TO  99  

ISKIP  =0 

ZEROl  = CMPLX  (0.,1.)  _ 

1 COSPHI  = COS(PHI) 

_30 SINPHI  = SIN  (PHI) 

^052  = COSPHI  • COSPHI 
SIN2  = SINPHI  * SINPHI 
OMFGA  r a » TOI 
FTA0=1. 000338 
35  AKP=OMFGA/C*ETAO 

OMFGAE  = OMEGA  • EO 
EPSIL=1.E-15 

10  AlO  = CMPLX(-C0SPHI,  0.) 

All  = Ain' 

40  Al?  = A17  = A20  = A34  = A35  = C mplX ( - 1 . , 0 . ) 

A56  = C«PL  X( 1 ., 0 .) 

A21  = CKPLX  (l.,n.) 

99  AK2=AKC*CSORT (CMPLX(E22, -SIGMA2/0“EGAE) ) 

ETA2  = AK2/AK0 
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45  ETfl22  = ETA2  • FTfl2 

TFMP  = ETft22  - SIN2 
7£Tfl?  = CSn^T  (TEMP) 

n44  = 7ETfl2  *(  SIN2  / (AMU2  * ETft22-SIH2)  » 1.)  /flMU2 

A45  = -A44 

50  AK3  = AKO  * CSORT  (CMPLX(E23,  -SI'^MAJ/OMEGAE)) 

ETA3  = AK3  / AKO 

FTA32  = FTA3  * FTA3  _ 

temp  = FTA32  - SIN2 

7ETA3  = rSOPT  (TEMP) 

55  A66  = 7ETA3  * (SIN?  / (AMU3  * ETA32  -SIN?)  + 1.)  / AMU3 

_ 70_CONTINliE  _ 

' AKl  = AKC  *CSQRT(  FMPLX  (F?l,  -SIGMAl/  OHEGAE)) 

20  ETAl  = AKl/AKO 

ETAl?  = <-TAl  * FTAl  _ 

60  TEMP  = FTA12  - SIN? 

7ETA1  r CSOPT  (TEMP)  _ 

A22  = A?3  = -7ETA1  *(SIN?/  (AMUl  * ETA12-SIN2  I ♦!.)  7AMU1 
A?3  = -A?3 

RCTAl  = AKO*  7FTA1  * 71*7ER01 
65  30  A3?  = r'^XP(-BETf.l) 

A33  = CEXP(BeTAl) 

A42  = A2?  * A3? 

A4'3  ='  A23  * A33‘ 

RFTAl  = AKO  * ZETA2  * 72  * 7ER01 

70  A54  = CEXP  (-RETA1) 

A55  = CEXP  (RETAl) 

464  = A44  * A54 

A65  = A45  * A55  _ 

T2'n=A54*A66-A64*A56  ' 

75  T 2n=A55*A66- A55* A56 

A1  = A2=A  3=A4=4RS(SI3MA2-SIGMA3) 

IF( ((Al.LT.FPSIL  ) .AND.  (42.LT.EPSIL  )).OR.  

1 ( (AT.LT.EPSIL  ) .AND.  (A4.Lf .EPSIL  ) ) ) GO  TO  31 
T 2=T?N/T?0 
8 0 (iO  TO  3? 

31  CONTINUE 
T 2=0. 

32  CONTINUE 

tTN  = ( A32  ♦ A44  - A4?  • A3'4)  - <A3?  * "A45  - A4‘2  • A3'5  ) * T? 
95  T3D  = (A33  * A44  - A43  * A34)  - (A33  • A45  - A43  * A35)  * T2 

A2  = A4  = A 0S  ( S IG  MA  ? -SI  (SMAT) 

A 1= A 3= A BS( SIGMAl -SIGMA?) 

IF( ((Al.LT.EPSIL  ) .AND. (A2.LT.EPSIL  )).OR. 

1( (A3.lt. FPSIL  ) .AND.(A4.LT.EPSIL  ) ) ) GO  TO  34 
go  IF(A1.\T.EPSIL)T3  = T2*A3‘2/A33 

IF(A1.GE.EPSIL)G0  TO  36 
R1=REAL (T?) 

B2  = AIMAG( T?) 

B3=REAL  (T3) 

95  B4=AIMAG(T3) 

115  FORMAK*  *,8  (E13.5,?X)) 

36  CONTINUE 

IF(  41.LT.EPS  IDGO  TO  35 

T3  = T3N  / T30 

100  GO  TO  35 

34  CONTINUE 

T3=0. 

35  CONTINUE  

RENUM  = (410  • 42?  - A?0  • 412)  - (AlO  • 423  -4  20  * 413)  * 13 
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10? REOEM  = (ail  * A??  - m * fll2>  - (flit  * A33  - ftZl  * ai3  ) * T3 


RE  = REMUM  / REOEM 
T£MP  = (1,  - RE)  / 

(1.  ♦ 

RE) 

PMI  = TEMP  * COSPHI 
50  AMPIMP  = CABS(PMI) 

110 

FAZIMP  = GANG  (PMI) 
IF  (FA7  IMP.LE.PI2I 

FAZIMP 

= FAZIMP  » PI 

IF  IFA7IMP.  GF.  .PI2»  FftZIMP  = FAZIMP  - PI 

LT.3.05  0)  FAZIMP=0.35  0?IF<FAZIMP.GT  .1.<4D0>  FAZI'^P  = 1.  400 
IF(  AHPIPP.LT  .0. 0 01»  AMFIMP=a  .001  JIF(AMPIMP.GT.  0.200  ) AMPi'1P=d,  200 


115  AMOIMP=1000.*AMPIMP$FA7INP=1000.*FAZIHP 

IFlFAZIMP.ir .0., OR. AMPIMP.lt. 3. IPRINT  1776 

1776’  FCRMAT(*0  AMP  OR  PHASE  NEGATIVE*) 

P REIU^ 

I FNO 


FUNCTION  GANG 


1 

r* 

FUNCTION  CANG(Z) 

c 

COMPLEX  7 

5 

c 

c 

c 

c.. 

...INITIALIZE  GANG  FOR  QUADRANT  CORRECTION,  GET  RE<7)  AND  IM<Z), 

TEST 

c* . 

POP  ♦ OR  - TO  FTNO  CORRFOT  HA L F- PL  ft NF - 

10 

DATA  (C  = 2.997925E8),  (T PI=6.  2831 8530 71 796) 

PI2  = TPI 

PI  = TPT  /?. 

PI02  = TPI  / 4 

CANG=0. 0 

16 

X=REAL(  Z) 

V=AIMAG(Z) 

c • • 

IF  (X)  10,50,90 

•••X  •LT*  O^Oc*************************************************** 

• • • • 

20 

10 

20 

IF  (Y)  20,30,40 

CANG=-PI 

30 

GO  TO  90 

CANG=*P I 

40 

RETURN 

CANG=+PI 

25 

c. . 

GO  TO  90 

50 

60 

IF  (Y)  60,70,80 

CANG=-PI02 

30 

C. 

RETURN 

. .X=0  AND  Y=0  IS  REALLY  UNOEFINEO,  BUT  IN  AGREEMENT  WITH  COC  HE. 

• • 

c. 

70 

. .RETURN  A VALUE  OF  CANOO.O 

CANG=0. 0 

• • 

60 

RETURN 

CANG=fPI02 

35 

c • • 

...y  -r,!.  n.n  and  y .ft.  n.n  ad.iu'^tft  for  rnPPFCT  nuAnpawT 

90 

CANG=ATAN  (Y/X)  4-  gang 

RFTURN 

END 
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SU3R0UTINf  SETUP 


r 


1 

SUBROUTINE  SETUP 

COMMON/ TB REE /F, PH  I, SIGMA  1 , E21 , A MUl , Z 1 ,S IGM A 2 , F2 2 , AMU  2 , Z2  , S IG M A3 , 

• 

E23,AMU3 

COMmON/COOE/SLCD(6  5, 3)  ,Ga  ,GLC0(14  8) 

5 

DATA  CFM/l.E-Z/ ,CSM /5  . / , CWG/ 1. E- 2/ ,C0G/1.  E-4 / 

TPI  = B.*  ATANd  .0  )IF=  l.E5$PHI  = 80.*l  TPI/360.1 

flHU 1=AMU?=AHU3=1. 

00  1 1=2,6'; 

1 

READ  1C  ,SLC0(I,1)  , SLC0(I,3)  ,SLC0(I,  2) 

10 

13 

FORMAT (10X,F10.1,F10.1,Fia.l» 

D02I=l,19fIB=(I-l)*8+l$IE=IR*7$IF(I.EQ. 19) IE=148 

2 

REACH,  (GLrO(K)  , K=IB,  IE) 

11 

FOPMAT(8F10. 1) 

DO  3 1=2,65 

15 

3 

PRINTS, I,SLC0(I , 1) ,SLCD (I, 2) 

SLCOd,  1)  = 1./SLC0(I,1) 

5 

F0RHAT( 15 ,2F15.1) 

PRINT6f0n4I=l, 148IPRINT5, I,GLCn( I) 

1* 

GLCO(I)  =1. /SLCOd) 

20 

6 

FORMAT(*0  GEOLOGY  RESISTIVITY*,//) 

GLCO(O) =SLC0(1,1) =CFWSSLCO<1,2)=20. 

SLCDI1,3) =20. 

Sir 0(65 ,l)=0FHJSLCO  (65,  2)=20. 

SLCn{65 ,3)=20, 

25 

END 

5.  DESCRIPTION  OF  PART  I 
5.1  Overview 

The  lithology  and  terrain  of  the  LORAN  coverage  area  are  digitized  from  maps 
in  PART  I as  illustrated  in  Figure  9.  The  required  three  maps  are  those  which 
contain  information  on  the  type  top  soil  or  overburden,  type  or  age  of  basement 
rock,  and  ground  elevation  above  sea  level.  The  first  type  of  information  is  ob- 
tained from  soil  maps  usually  published  by  the  country  of  the  interested  area. 

These  are  colored  maps  with  a legend  defining  the  various  types  of  soil  and  associ- 
ated depth.  Geological  maps  are  available  for  most  areas  of  the  world  and  the 
colored  legend  indicates  the  type  of  rock.  Terrain  maps  with  elevation  data  are 
readily  available  for  the  entire  world.  Fortunately,  these  have  been  previously 
digitized  by  the  Defense  Mapping  Agency,  and  magnetic  tapes  with  this  data  in 
proper  sequence  are  available  upon  official  request. 

The  following  equipment  was  required  for  data  digitization: 

Calculator  - HP9830A 
Digitizer  - HP9864A 
Cassette  Memory  - HP9865A 
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■ 1 


AND  DEPTH 


j Figure  9.  Formation  of  Data  Base  - Part  I 

I 

I 


The  maps  are  digitized  by  recording,  on  a cassette  tape,  the  contour  table  coord- 
inates in  inches  using  the  southwest  corner  as  the  origin.  The  cursor  cross  hairs 
are  moved  around  each  constant  color  contour  and  points  are  generated  in  incre- 
ments of  one  hundred  to  the  inch.  The  table  coordinates  and  geographic  coordinates 
of  the  four  corner  points  of  the  map  are  recorded  for  later  use  in  conversion  of  the 
former  to  the  latter.  These  cassette  tapes  are  transcribed  onto  a disc  file  by  way 
of  the  CDC  6G00.  The  merged  cassette  data  is  then  copied  on  to  a 1/2-in.  2400-ft 
800-bpi  tape  reel.  The  data  on  this  tape  is  converted  to  LAT/LON  and  color  num- 
ber, and  then  into  scan  lines  by  PROGRAM  GGCI.  The  scan  line  output  format 
contains  24  bits  for  longitude,  24  bits  for  latitude,  and  12  bits  for  color,  making  up 
the  60  bit  CDC  word.  These  points  will  then  be  sorted  in  PROGRAM  SOR  with  the 
primary  key  being  the  latitude,  secondary  key  the  longitude,  and  tertiary  key  the 
color.  PROGRAM  DECODE  sequences  the  sorted  line  points  into  a 30  arc  second 
matrix  lor  the  final  output  tape  of  PART  I.  These  programs  ar“  listed  in  Section  B. 
Three  such  magnetic  tapes  are  required,  one  for  each  of  the  ground  properties. 

For  a coverage  area  of  48  square  degrees,  there  will  be  almost  700,  000  data  points 
recorded  on  each  7 -track  binary  magnetic  tape. 

The  Defense  Mapping  Agency  generates  similar  data  tapes  in  a program  titled 
"I.ineal  Data  to  Terrain  Matrix."  This  program  is  capable  of  forming  a data  ma- 
trix from  map  contour  scans  with  a three  arc  second  to  a one  minute  spacing.  By 
using  the  standard  DMA  format  as  input  specification  to  PART  II  of  the  prediction 
program,  DMA  and  locally  generated  tapes  can  be  used  interchangeably.  The 
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present  program  is  setup  for  a thirty  arc  second  grid.  Changes  in  this  spacing 
will  require  changes  in  the  incremental  sequencing  and  storage  allocations. 


3.2  l‘rut;ram  l.istin);  for  Part  I 


PR0r,(?#M  SOR  (TAPE3,  tape?,!  APE*. , TAPEl  I 
9 READ(l)  IP 


— iFteorrm?,  I 

1 HRITEt  2) IPJGO  TOP 


CALL  SMENO 
REMIND  3 


— 2096 


•*4-'COHISH 


9— RFAntrrTP 

IFtEOF (3)»4, 3 

--  

4 STOPJEND 

*4C60,T160,CM77777,IP2.  

FTN,SL,A,P,eL,R=3,FL=20300. 

VSN,TAPE1  = 0S  0139. 

RF  GUEST,  TAPrr,~MT,f,ITO. 

WSN,TAP£2=OS017?. 

FEaU£ST,TAPE2,HY.L,RI-N£l^  _ . 

FILE  (TAPE?, RT=S,0T=C) 

ldset(files=  tape?) 

L60-*  - . - -- — 

PROG PAM  OECOOEt INPUT, CLT PUT, TAPE 6=0UTPUT, TAPE 5=1 NPUT, TAPEl, T APF?) 
COMMON  I0UT( 721)  ,ICURPT, IM INH, I M AXH , IM INH2 , IMAXW2,INFL00 

DIH»-NSI0N  T(  9001!  , ir.  (5C01)  ,TpUF(912) 

OIMFNSICK  M0NE<7  2i),MT).0t721),MSAVE(721) 


11 


ISTO'P=  0 
XXOE3  =6. 

DO  11  I6=l,7?l 

lOUT (I  6) =0 


"r?o=5r3- 
LEN=51 ? 
1T(3T  = 0 

XCUR=0  . 

xc»a-x-xoE& 

IMINH=  0 
IMA  XM=  0 


IHTNW? -"lOOB 
IMAXH2=-1000 

- 

IE?=0 

C-SKIP-  IES^'L-OA-TA— 
IPMASE=1 
15  CONTINUE 


IF(IF?.F0.1)  TPHASE=? 
J=1 


^IPST  = 1 


CONTINUE 
INFLOO  = 0 

4 TS  TO  READ  IN  ANOTHER  SCAN  L TWF 


IF( IFIRST.EO.l)  G0T06 
IF(XCUR.LT. (XXOEG-,004))  CD  TO  6 

•*  -OECOOE  MOOIFICA-IIONS  ....  .(SIARI) 

**  -MODIFICATIONS  (BRUCE  HAINES  8/18/76) 
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00  aiToi>  KA=i  ,r?i 

a005  MSaVE{KA>=IOUT(KA) 

IF( MCA SE.GT. 1)  GO  TO  8010 
no  BDD7'KA=1,72I  ‘ 

8007  MONE(KA)=IOUt(KA) 

CO  TO  ^ 

8010  IF(MCASE.GT. 2)  GO  TO  8030 
no  8015  KA=1,721 
8015  MTWO (< A)=IOUT (K A ) 

GO  TO  8250 

c~  ■ **  -rcR  NOpysL  ?cans 

8030  00  8100  KAil,721 

IF «MTMO{<A). to. 0 ) GO  TO  8C50 
I0UT(KA)=HTW0(KA> 


— f,e  TO  8100 

8050  IFdOJT  (KA),  EO.O  ) GO  TC  80  70 
GO  TO  8100 

8070  TFrFTONrtKAl.EO.P)  GO  TO  3U80 
IOUT(KA)=MOME(KA) 

GO  TO  8108  - - 

8080  IF(KA.EQ.l)  GO  TC  8100 
IOUT(KA)=IOUT(KA-l) 


8108  COWTINOe 

00  8105  KA=1,721 
MONE(KA)=HTWC(KA) 

HTWOfKAT^SAVTTKAr 

8105  CONTINUE 

GO  TO  .8388  


**  -FOR  LAST  SCAN  LINE  ONLY 


8170  00  820  0 KA=1 ,721 

IF <MTWO<KA) . EO. 0)  GO  TP  8180 

rour(KA)  = MTwc(Kin 

GO  TO  8200 

»488-3F  tMOHElXAl.EQ.Q)  GO  30-8108 
lOUT (K A)=MONE (KA ) 

GO  TO  8200 


— 8150  TFTKA.eO.l)  CO  TO  8280 
IOUT(KA)=IOUT(KA-l) 
,8200  


GO  TO  8300 


**  -FOR  FIRST  SCAN  LINE  ONLY 


8250  00  8270  KA=1,721 


UU 


82  55 

.JNE(KA) 

GO  TO  8270 
8255-  IF  tHT<<  or  KAl  _€0-81-  00^-10  -8200 
lOUT (KA)=MTWO(KA) 

GO  TO  8270 


8260  TFtKA.EQ.H  GO  TO  82T0 
IOUT(KAl=IOUT(KA-l» 

^8273  continue  „ 


8300  HRITE(6|612>  XCURS  

612  FORMATax  ,G21.1A) 

WRITE(2)  XCURS 

WRITgt  6.610)  TOOT 

610  FORMAT  (2514) 

WRITE(2)  ( IOUT(  KLH)  ,^KLM=li721) 

~8i.Trff  5n:uRS=  xcirR  ■■  ^ 

MCASE=MCASEe1 

1F( TST QC . EQ . 1>  STOP  - 

IF (XCUH. GT. ( XXDEG+.991) ) GO  TO  8170 


6 CONTINUE 

XCNT=K  CNT+ (1 ./130. ) 

tXCKT.LI.<XCUS*<  ) GOTO  4 

IF(IP^ASE.EO.^)  GOTOS 
IF(IE? .EQ.l)  GOTO  15 

-9 If  HE.gP.l)  eOTO  9» 

no  9 15=1,5001 
Y(I5)=0, 

IC(T5»  =C 
9 CONTINUE 

IFdFIBSI.eO.ll  &0T07 
XCUR=XNEW 
Y (1>=YSAVE 

TC  m-Te?ATg 

7 CONTINUE 
20  CONTINUE 

no  10  rT=j,  50  01 

IF( 13. LT.500 1)  GOTO? 

Pf?INIGC3»XCUfi,  YTSOJD,  10(5000) 

603  FORMAK*  ARRAYS  Y ANO  IC  TOO  S N ALL  • , ?&?1  . 1 4 , 15  ) 
STOP 

2 CONTINUE 

if(ip4ase.ne.i)  goto  •'1 

RE  An(5  ,*)  XNE  W,  Y <13  ) , in  (I  3) 

IF(FoF(5))  49,3? 

C READ  IN  A POINT, 

31  CONTINUE 

RFAO(l)  IP 


1 XNEM=FL0AT(SHIFT  (IP,-36))/120, 

TF (XnER.LT.XCURI  G0T031 

Y( I3)=FLOAT<SHIFT(IP.-12) .ANO.  7777  7777  R)*.001 
IC<I3»  = IP.  ANO.  7777  6 
IF(IC( 13)  .GT  .300)  GOTO  220 
C FOR  GLITCHES 

TFTTXNEW.tT. 9).0R.  ■ - ■ :T 

IF((Y(  13)  .LT  .47)  .OR.  ( Y ( I 3)  . GT  . 5 5 ) ) GOTO  221 

_32 iF(J.EQ,l,AND.i3,E0.1)  XCL'RrXNEW 

TF(XCJP. NS.XKFwr  GOTO  J 
IF(I3.E0.1)  GOTO  10 

IF  TYTUT  .4.1,  Y(IJ-lTi  OOTO  302  - 

10  CONTINUE 

PRINT  e04 


STOP 

-25  P^^^T*  STMT  NO.  96  ENCOUNTEREO.  FINTSHFD.” 

220  PRINT*,"  IC.GT.300,  PROBABLE  GLITCHES." 

GOTO  31  

221  WRITE! 6,2601 ) 

2601  FORMAT!*  COORDINATES  CUT  OF  RANGE.  PROBABLE  GLITCHES.") 


< u ox 

302  WRITE! 6,3601) 

3601  format (6 (-  A point  out  of  ORDER"! ) - 

Y( 13)=Y< 13-1) 

GOTO  10 


50  continue 

P^NT*,"ENO  OF  OATA.  ITOT=".ITOT 

C IE=ENO  OF  FILE  SWITCH 

GOTOl  


IFIRST=0 

^ ,L=13-1  . ^ 

C GO  From  1 To  L with  ARRAYS 
YSA  VE=  Y(  13) 

ICSAVE  = ir:<I3) 

Y ( I 3 ) = 0 . 

IC( 13) =0 


IF(XCJR.LT. ( XXDEG-. 004) ) GO  TO  4 
C PRINT  scan  line 
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1 r i = i x t i u 

INFLOO=0 

??fU^(^(ici}fiPT«-i>-y(iruoFT)).i.E..oi';i)  g oxer  51 

C IF  WITHIN  .015,  GO  TO  51 
GOTO  76 

51  IF(IC( ICURPTtl) .EO. ICURCL)  G0T077 

ISP1=0 

GOTO  TT 

77  CONTINUE 
^gOPPT^ICURPT.l 

76  IF( ICURPT.GE.L)  GOTO  4 

73  ICU«PT=ICUOPT*l 

IF( IC( ICURPT ) .NE  .ICUOCL)  GOTO  45 

CALL  F ILL (Y ( ICURPT-l) , Y ( IC LRPT) , IC ( ICUPPT) ) 

IF( APS (Y ( ICURPT+ 11 -Y(I CUPP T)).GT.. 0151)  GOT074 
C IF  NOT  HITPIM  .015,  GO  TO  74 
ICURPT=ICURPT+i 

IFdCdCURPT-ir.EQ.ICdrUFFT))  GOTO  75  

102  ICURCL=IC (ICURPT ) 

GOTO  2 5 

75  CONTINUE 
GOTO  102 

— e 

45  CONTINUE 

ISC1=IC( ICURPT) 

IFTAP; (Y (ICURPT* 1)-Y (ICURPT)! .GT..D151)  GOT07B 
C IF  NOT  WITHIN  .015,  GO  TO  78 
GOTO  79 

78  IFdCURPT.GE  .L)  GOTO  4 
ICURPT=ICURPT-1 

GOTO  ru 

79  CONTINUE 

82  ICURPT=ICURP T+1  

" IF(TCU17CL.NF.IC'drURPT7T^(lTmnrD 

G0Tn81 

- 40  CONTINUE  

IF( < I99H1.EO .I99H2) .AND. (T99H2.EQ. ICUPPT) ) GOTO  105 
I99H1= I99H2 

I99H?--ie»RPT 

106  IF  (ISP1.N9. 0)  GOTO  83 
ICURPT  =ICU9PT-1 
ICURCL'=TC  (ICURPT! 

GOT  025 

—105-  ISP4=ISP4*1 
GOTO  1 06 

83  ICUPPT=ISP1 

ICURCL  = TSC2 

G0T075 

51  ^ a|j^^F^}^^  Y atURPT-ai  iTlITiLREl!  iIC_<  ICUPETJJ 

G0T025 

C ~ 

72  CONTINUE 

e MOLC  sssuitoTCAsri) 

74  CONTINUE 
'C  IS 

C 4 IS  TO  READ  IN  ANOTHER  SCAN  LINE 

ICU8PI=-ICURPT*1  

ISP1=ICUPPT-1 

CALL  FILL(Y(  ICURPT-l)  ,Y(ICURPT)  ,IC(ICUPPT)) 

rrtAqS-tTTTCURPTYlTT  (TCURPT)  1 .GT.  .0151)  GOTOT? 

C IF  NOT  WITHIN  .015,  GO  TO  72 
ICUPPt=ICyRPT*l 
tSC?=IC(ICURPTt 
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C (IN  C4  SE  B flCKHARDS) 

ICURCL-ICIICURPTJ  

GOTO  ?5 
ENO 

SUBBOUTINf  fILL — (YNT  N , YWfty-,Tei 

C FILLING  IN  COLOR  NUMBERS 

COMMON  I0UT(7?1) jICURPT, IM JNH.IMaXHi IMINH2, IMAXH2 , INFLOQ^ 
INFLOO  =INFL5o+1 
IF  ( YM  IN, LT. 48.  ) YMIN  = 4fl. 

IftYHAX  .Cl,  54,)  YMAX=^54. 

IF ( YMA X. LT.48. ) YMAX=4e. 

IF(YMIN.GT.54.)YMTN=54. 

IF  (YMTN. ST. YMAT)  GOTOl 

IARYMIN=(YMIN-48.>*1?0.+1.+.5 
lARYMA  X= (YMA  X-4  8 .)*l|0,  + l,+,5 
C IT  IS  ASSUMED  T4AT  THE  REAL  VALUE  HILL  BE  TRUNCATED 
IF(( IARYMIN.EQ. IMINH) . AMO. (IARYMAX.E3. IMAXH) 

- 2 - ^Nn,(HaN*i,£Q,IMI-NHE)  .ANO.  ( IMAXH.E  0.  IMAXH2))  GOTO  AO 
31  IMINH2=IMINH 

IMAXH2=IMAXH 

IMINM- lAB YMIN 

IMAXH=  lARYMA  X 

C FILL  UP  AN  ARRAY  OF  LENGTH  721 
■"  DO  1001  I=IARYfnTr,TARYHAX 

lOUT (I )=IC 
CONTINUE 

return 

1 PRINT51 


-9t- 


FOBMAr(*  ABORTED  MECAUSE  SUBROUTINE  HAS  CALLE-B-Hf-T-H  YMINtYMAY*) 


STOP 

30  ICURPT =ICURPT+3 
Ga^TD  'TlSENn  - 
MC64, T450, CM  25500 0 ,TP2. 
-FTNt  $4^4=99  990) 
VSN,TAPE1=GGC64C. 
REQUEST, TAPEl.HY.L. 
VSN,TAPE2=OS0139. 
REQUEST,  TAPE2,HY,L,RING. 
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MCCOMISH 


m. 


LI)T,CC(1CI), 


EXIT.S. 

REHINOtOUTPUT.  — 

REWIND, TAPES . 

COPYCF, TAPES , OUTPUT. 

PROGRAM  SGC 1( INPUT, OUTPUT .TAPE 8, tape 1 ,T  APE?) 

COMMON/INPUT/IRECNO.SUMVEC , I XS T , I Y ST  , t V FC ( 6 Q 0 0 ) , IRECTYP, IFEATNO 

CQNMON  IBLKNOiISTVEC.IAGE. lNOLO(i06P) .jRECRES.LiNOFTR^C 

COMMON  XP.YP.RCP  (5  0)  ,SlO^R(20)  ,A(10),P(^10),AA(10)  ,8BTl( 

7 00(10  ) ,04850,05052.05254 
DIMENSION  NOXLAT(6000) ,NOVLON&(600C) ,LAT (6000) ,LONG(&JflO> 
OIMENSION  NXX  (5  ) 

DIMENSION  R(6000), 01(600 O.NGEO (6000), NO VEC(150), 02 (6000) 

DIMENSTOW  TCCNT(8030)  , SL  OF  EtStTO  OT 

1 FORMAT  (14) 

2 FORMAT (•  NUMREC=*I4) 

^ 94  F0RHAT(  5(3X»L0nGIT00E  LATITUDE*3X) ) ~ 

106  FORMAK*  TAPE  HEADER*) 

_ RECORD  NO,*  t4,6X-,t66ATU«E  NO. *I4,6X, ♦RECORD  TYPE»I4.6X.' 

ARLOCK  N0.»I4,4X, *RECOROING  RESOLUTION=*  13) 

116  FORMAK*  RECORD  NO. *14, 6X,*FEATURE  NO  . * 14 , 6X  , * RECORD  T YPE  * 14 , 6 X , ' 

IFLOCF  N0-.*T4) 

135  FORMAT(*  AGE=*I12,4X,  13*  RECORDS  IN  THIS  FEATURE*) 

145  .format  (*  first  X=*f6.gX. , ^._^*FIRST  V = *16,2X, 

• niim  uF0  = ^AT7y  .»^STARTTnG  v#FC.  Pns.  N0.=*I3) 


NUMVTD^*T572X,»^STARTTnG  VEC. 

, ERROR  IN  VECTOR*  15) 

F0RHAK^R44X,E12.^T)-  - - 
921  F0RMAT(7(2X,F11.3)) 

930  FORMAT (7 (2X, 16, IX, 1 6, IX,  13)) 


371  FORMAT!* 

— 


-974 FORMAT  (*- 


935 

93£l_ 


lNDPOTNTS  uf  vfctdrs 


FORMAT (/*  TOTAL  NUMBER  OF  POINTS  =*I5) 
_EQRMAJLtS_l5X,_*LQNG*6X*  L AT*7  X)  1 


T=rONGTTUDEi 


LAT.»T 


938  FORMAT (12 (2X ,13 , IX , 15) ) 

944  FORMAT (♦  NUMBER  OF  POINTS  IN  FEATURE  MAY  EXCEED  ARRAY  SIZE*) 
950 — FORM  AT  (4(  1X,E8.  C-y  1XtF6  ,6, 1 X-,  F 7,0 , 1 X , F 7.0)) 

951  FORMAT  (10 (IX  ,15  , IX  ,16)  ) 

952  FORMAT ( 6(2X,  020)  ) 
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9W 

PHIR1=48.0  £ PMIR6=50.0 

THETA1=6.[1  £ TH£TA2  = 8.0  £ 

COWPOTATTON  ♦) 

£ PHIR11=52.0  £ 

TH£TA3=1Q,  £ THETA4=1Z, 

PHIRl 6=54. 0 
£ THETA5=14. 

READ  1,N(IMPEC 

PRINT  B.NUMREC 

Mnf  11  E -A 

IRECNO=  -1 

DO  60  8=1,150 

eir 

70 

NUVti;)  M)='J 

CO  80  K=l,6000 

NDXLArK)  = Q 

80  NOYLCl'jC(j<)  =0 
C READ  IN  NEW  RECORD 
100  CALL  INGGCAI IEOF,NOWORO) 
IF(IEOF.EQ.O) 104,102 
10  2 NOFILE=NOFTLEM 

104  * 

105 


IFdRECNO.GT.O)  107,105 


ills 

107  IF (IRECNO .EQ  .1) 108, 110 
100  SLOMF  G = SLOWER!  5»  4 

SLOME  7=SL0MER(  H S 

SLOME  8=SLOWER<  6)  £ 

STOWE g=SLnWEPt  21 1- 

SLOME  10=SLOMER(  7»  £ 

SLOME  ll=SLOMER(  3)  $ 

SLOMP  12=SL0WFR(  8)  £ 

SLOME  13=SLOMER(  4)  J 

slowf  i4=slomer<  q>  i 


AG  =A< 

A 7=A ( 1) 
A 8=A(  6) 
A q-At  2) 
A10=A(  7) 
A11=A(  3) 
A12=A(  8) 
A13=A(  4) 
A 14= A! 


4 

£ 

£ 

“£~ 


q> 


5 G=D(5> 
P 7=6(1' 

6 8=B(5) 

B10=B( 7) 
611=6(3) 
B12=B(3) 
B13=P(4) 
BlG=8<q> 


110 

111 


120 

122 

130 


140 

141 

142 


RCP4  8= (RCP< 1) tRCPC  2)  +RCP( 3) ♦PCD (4) ) 74. 

RCP5  0=  (RCP (5 ) ♦RCP(G) ♦RCP (7) fRCD (8) ♦RCP (9 ) )75. 

PCP82=  tRCPtiOiY-RCPtm  *pcrt-r?)  ♦rcphs)  ♦rcp(ib))/8. 

BCD  84=  (RCP ( 1 5) ♦ RCP (16) ♦ RCP ( 1 7 ) ♦ RC P ( 1 8 ) fRC P ( 1 9 ) ) 75 . 

^PRINT  jj,^i^OM£g^S^OME7,SLOMF8,S-  OME 9 , S LO ME  1 0 , SL OMEl  1 , S LOMEl 2 , 

PRINT  921  , (RCP4  8 .RCPFO  ,RrF5? ,POP5  4,0485  0,05  052,D5254) 
IF<IRFCN0.t0.2) 111, 114 

PRINT  112,  IRECNC,IFEA  TNO,IREC  c'YP  ,I  9LKNO  ,IRECP.-:S 
GO  TO  120 

IPECTYP-,  ■ TPLKNO 


118  , 1 cuniu  , 

IFCIRECNO.GT .4) 122,100 
IF(IR£CTYP,EC,30)130,140 
NT=0 

PRINT  135,  IAGE,NOFTREC 
lEtNOFIREt.GT.G) PRINT  944 
GO  TO  100 

IF  (I9L  KNO.EO .1) 141,  142 
NOYLAT(l)=  IXST  £ NOYLONG(1)=  lYST 


PRINT  145,  IXST, 
N = NUM(/EC 

DO  600  1=1, N 
J=  I ♦NT 
K=ISTVEC  ♦ I 
IYFC( J) =INOLn(K) 


lYST, 


NUMVEC, I6T7EC 


GO TO (2 10, 220, 230, 240, 250, 260, 27 0,280, 2 90, 300,310, 320, 33 0,340, 350, 3 
660 ,370  ) JVEC 

C DELTA  X IS  THE  INCREMENT  OF  LATITUDF  FOR  EACH  SCALED  VECTOR 
C DELTA  Y IS  THE  INCREMENT  OF  LONGITUDP  FOP  EACH  SCALED  VECTOR 
210  inELX=0  $ IOELY=2  £ GO  TO  500 

220  IOELX=l  £ I0ELY=2  ? GO  TO  500 

“ * inFLY=2  £ GO  TO  500 

TOEtY=l-^t;tr  T-Q-WDO 

IOFLY=0  I GO  TO  500 

IDELY=-l|  GO  TQ  5C0 


230 

-?tr(i- 


250 
260 
27  0 
280 
290 
300 
310 
-rttr 
330 
340 
350 
360 
370 
372 
500 


IOELX=2  £ 
IPCLX=? - 
I0ELX=2 
IDELX=2 
I0ELX=2 


TDELY=-2f  GO  TO  500 

I0ELV=-2?  GO  TO  FCO 

IOFLV=2  ? GO  TO  500 

_ _ ICELY=2  £ GO  TO  500 

IDELX=-2$  IDELY=1  f GO  TO  500 

“ I GO  TO  50  0 

IOELY=-lf  ■■  - “ ■ 

I0ELY=-25 
IDELY=-2? 

IOELX=0  $ I0ELY=-2£  GO  TO  5C0 
PRINT  371,4 

IDELX=0  £ IDELY=0  £ GO  TO  5C0 
NOXLAT  (J^l)  = NOXLAT  (J)  ♦ lOELX 
WOTLOiIGTJtI)  =NOTLONGt  J)  ♦ TOFLY 


IOELX=l  ^ 
IOELX=-lj 
IDELX=-2( 


I0ELX=-2S 
IOELX=-25 
IDE  LX  =-l 


GO  TO  5Cn 
GO  TO  500 
GO  TO  500 


67 


^600  CONTINUE 

NT=NT+NUMVfC 

IF ( IBLKNO.lt .NOFTPFC) 100,602 

682  NOW£C< IA6€>=N0W£C<IA6F»*  WT  -- 

NT2  = NT  +1 
PRINT  CJL 


PRINT  930 , (NOYLCNGC J) , NDXL AT (J) , lAGE ,J =1 ,MT2 ) 

, PRINT  936, NT2 

^ DO  701  J=1,NT2 

: ICONT  IS  THE  PACKED  ROPD  OF  COORDINATES  IN  X-Y  UNITS!  LONG. =29  BITS,  L 

: 29  BITS,  AGE=12  BITS,  ALL  RIGHT  JUSTIFIED 


R(J)=SCRT(  ( NDXLAT(J)  -XP)*»2  + ( NOYL ON G ( J ) -YP ) • * 2 ) 

IF(P.{J).GI.RCP98>LAT(J)={98.-2,Q*(P(J1-RCP98)/(D9850  J)*lpp0, 

IF(P(J) .GT.RCP50 .And. R ( J ) . LE . Pfi P98 ) L AT ( J ) = ( 98 . 0 > 2 . 0* (RCP98-R ( J> ) / 

A (09350  ) )*100C. 

IF<R< J >,GT,RCP52.AN0,S( J» ,L£.RC050>LAT  < JT-< GA. A*2,6*  <RCP50-R ( J) ) 7 
A (05052  ) )»1000. 

IF ( R( J ) , GT.RCP59.AND. R( J) . LE . RCP52 ) LA T ( J) = ( 5 2 . 0 +2 . 0 » ( RC P52 -R ( J ) ) / 

A » -)•*  ll?0-(h 

IF(R(J  ) . LE.RCF59  )LAT  ( J)  = (59  .+2.  0»  ( RCP  5 9-R  ( J)  )/'(D5259  ))»1000. 

SLOPE! J) = (YP-NOYLONG(J)  )/  (XP-NOXLA T ( J ) ) 

IF(SL0PE(J) .GE. SLONE  6 .AND.SLOPE(J).LT.SLOHE  7)  611,615 
Z different  FORM  OF  THE  POINT-SLOPE  DISTANCE  COMPUTATION  FORMULA--SANE  R( 
fell  OKJMABS!  (SLOME  6 ♦NOXtJlT(J»  - NOVLCNGIJ)  *A  fe  >/B  fe  ) 

D2(J)=ABS(  (SLOME  7 »NPXLAT(J)  - NOYLONG(J)  +A  7 ) /B  7 ) 

LONG(J)=  ( 6.  O + DK  J)  / (01  (J  I FD2  ( J)  ) ) *10  00 


615  IF(SLOPE(J)  .G£. SLOME  7 .A  NO  . SL  OP  E ( J ) .LT  .SL  ON  E 8 ) 

LONG(J)=  ( 7.0f01(J)/(01(J)+02( J)) ) *1000 
GO  ITT  700  - - - 

620  IF(SLOFE  (J)  . GE.  SLOME  8 . A NO . SL  OPE  ( J ) . L T .SLOM  E 9 ) 

621  D1(J)=ABS(  (SLOME  8 »NOXLAT(J)  - NDYLONG(JI 


LONG(J)=  ( 8 .OtOKJ)/ (Cl  (J»+02(  J)  > ) *1000 

625  IF(SLOPE(J» .GE.SLOHE  9 . AND. SL OPE (J1 . LT. SLOME  Iff) 

626  01(J)=A9S(  (SLOME  9 *NDXLAT(J)  . - NDYLONGCJ) 

02(J»  = ABS(  (SLOME  10  TNOXIAKJ)  - NOYLONG(J) 

LONG(J)=  ( 9 .0  + ni(  J) /(Ol  (J)  ♦02(  J) ) ) *10  00 

GO  TO  700 


616,620 
FA  7 ) /»  7 ) 
♦A  8 )/B  8 ) 


621,625 
+A  8 )/0  8 ) 


626,630 
♦A  9 )/9  9 ) 
♦A  101/B  10) 


02(J)=AeS(  (SLOME  19  *NOXLAT(J) 

^gNG^J  ) = P ( 13  .0*01!  J)  / (OK  J)  *0  2!  J)  ) ) ' 

65  0 IF  (SLOPE  ( J)  . GE.  SLOME  19  >651 ,655 
>51  01<J)=ABS(  (SLOME  13  *NOXLiT(J) 

D2(J)=ABS(  (SIOME  19  *NOXLAT(J) 

LONG  ( J)  = (19.0FO2(J)/(Ol(J)-n2(J)))' 

GO  TO  700  

655  IF(SLOPF( J) .LT.SL0ME6) 660, 693 

LONG(J)=  ( 6.0-0l(  J) /(02  (J)-Ol  (J)  ) )< 
60  TO  700 
693  PRINT  999 


- NOYLOM6(J) 

- NOYLONG(J) 
1000 


‘1000 


*A  13)/B  13) 
fA  19)/9  19) 


I 


fi 


VllB  0 0«li 


L0NC(J)=L0NG (J-1) 

700  NGEorjr=SffIFT(L6NG{  J)  ,36r.M.  SHlFTTLliT  tJ)  ,12)  .“OR.IAGE 

701  CONTINUE 

PRINT  g50,(R (J) ,SL0PE<J),D1(J),02{J) ,J=1,NT2) 


PRINT  951 , (LONGTJ) ,LftT (J), J=1,NT2) 

2^^  ^223^fffli?Ml)’ LffNlTT  ir,LaT(iT,  IXGf,NT2 

NXXC 1) =IFEftTNOSNXX( 2) =LONG(l) |N  XX ( 3 ) =L AT  < 1 ) INXX (4)  = I AGE  SNX X (5 > = NT2 
- 8UE:FgRWT<2rl)T>UXT14~rN*X4WirE<U»HT4^7^»1791,179W47^1 


buffer 0UT(2| 1} (NXXfll fNXX{5) ) fIF(UNIT{2i  11791 -1791i  1791 

17  91  BUFFER0UT<2f 1)  ( NGE0(  1)  ,NGEO( N T2 ») I IF (UNI T ( 2) ) 1792, 179 2, 1792 
1792  CONTINUE 


1000  CONTINUE 

^ _PR^l_aM.  (I  .NfiYECd)  .1=1. 1481 

SUBROUTINE  I NGGC 6 ( lEOF , NOHORO) 

THIS  VERSION  OF  INCCC  IS-FOR  THE  6 DfcGREE  TAPE 

COMMON/ input /IRECNO.NUMVEC , I XST, I YST , I VEC ( 6 0 0 0) , IRECT VP , I FE ATNO 
COMMON  IBLKNO, ISTVEC, lAGE. INOLD<1080) , IREC RES ,L , NOF TREC 


Z 00(10  ) ,04656,05052. C5254 


ITPt  REAL  LX,LT 
5 FORMAT (4(2X, P1E16.8I ) 

fe F0RMAT^<5  ( 2X  , PIF 1 T 

8 FORMAK*  OISTANCES  BETWEEN  PARALLELS*/5  (2X  ,E14. 8»  ) 

9 FORMAT (4(2X,P1E10. 4)1 


144  FORMAT (5(2X, 12, 2(2X,F7.1)) ) 


2 

436  FORMAT  (* 
—4^55^  FORMATI-*- 
496  FORMAT  (* 
1000  FORMAT (* 


ItlH.  t ,<:A,  > r 

FILE  SUMMARY  RECORD  *) 

4nvAlI0  record  TVP£s 

END  OF  FILE*) 

INVALID  AGE  NO.*) 


110  PRINT  1500 
120  IEOF=0 


IRECNO  = IRECNO*  1 



IFEATNO  = IN0L0(7)*16  ♦ 1N0L0(8)*  1N0LD(6)*256 

IRECTVP  = IN<}LOI9)  * (2**^)  * INOL0(40) 

IBLtCN0=  IN0L0(11I*(2**4)  ♦ INOLOdZ) 

IF ( IRECNO. EQ.l) 135,150 


PRINT  143 

IMT— 1 4-4.4 1 .4 X 1 1»- j4  Yd)  ,1^1.  t9) 

XPN={LX(5)*LX(19) )* (LY(15) -L Y (9) ) ♦ (LX (5) *LX (9) ) ♦ (L Y( 19) -L Y( 15) ) ♦ 

A (LX(9) *LX(15) ) *(LV(5) -LY(19) ) ♦(LX(15) *LX (19) ) *(LY(9) -LY (5) ) 
XP0=(LYtl5)-LY(5>)*(LX(19) -LX (9) ) - ( LX ( 15 ) -LX (5 ) ) • (L Y ( 19) -L Y ( 9) ) 
XP=XPN/XPO 

15)  » • (1  Y f 1 5)  M Y »<;)  ) / «LX  « 1 FI -1  X f 


Ml] 


i«xid 


1441  0L(M)=SQRT((LX(ME5)-LX(M)) **2  * (LY ( M+5) -LY ( M) ) **2  ) 


05  254=  <OL (10) *04 (11) *04 ( 12) *04 ( 1 3) *04(14) ) /5. 
PRINT  8,fOL(M),M=5,14) 


riAi.r,icuu,  AT,  If-  

DO  145  1=1,19 

145  RCP(i)=SQRT(  (4 X (I ) -XP) **2  ♦ (L Y( I) -YP) **2  ) 

PRTfrr  I46,d,RCP(I)  ,T=Tri9T — 

00  144  3 M=1.4 

Dt^GR  tMlaSAia<4(.X(MA5)-.TJ(aM*»))**2^  <LY  (W»5)  .LY  (H*( 


1<*43  PHI<‘I>  = ( ftSIN(DCP50  (M)/RCP  (H+4)  )»•(  180  ./3. 14159) 
PRrNT9,(PHI(  M)  ,t1=l,4) 

00  14T  H-t,9 

ft#  (M)  = YP-L V(M) 


i 9R(M)=LX(M) -XP  t CC(M)=XP*LY (4) -LX (M) »YF 


D0tM)=SQRT (flft(M) ’’a  ♦9B(M)»*2) 

= rn.x(K  )»Ypr-fLY{K  )*xp))/(lx(m 


)-XP) 
) -XO) )•»  2)  + 1.  ) 


SCMT  =cn.X(K  )»YPr-fLY{K 

147  BIM)=S0RT((( (LY (M  )-YP)/(LX(M 
DO  148  M=l,iq 

148  SLOMER (M)=(LY(M)-YP)/<LX(M)-XP) 

PRINT  5, (SLOMERl M) ,M=1,4) 

PRIwr-e,  tStOWgRfN)  ,11  = 8 ,19) 

150  IF(IRFCNO.EQ.2) IRECRES=IN0L0(23)*(2**4)  ♦ IN0L0(24) 

16  0 

162  JH  = IN0LD(45)  •(  2**4)  ♦ INOLO(46) 

JT-^  I K0L&<47»  » < 2**44  *--IWOU)X44*  - 
JU  = INOLO(  49)  *(2»*4)  ♦ INCLD(50) 

NOFTREC=IBLKNO  - 1 

IF  UU.  GT.  1T5  .aNCr.'.W.LT  .186)  170,  1T6 

^170  IF( JT. GT.175 .AND. JT.LT  .186)  172 , 186 

^ ^ IF(IA6E,GT.g.ftND.IflGE.LT|l4g)  50  0 tl8^ 

t74-  iFl  32^0«*gU»^C»4i4)^^-48i- 

177  IF(JT. GT.175.AND.JT.LT.186)17d, 180 

178  IF ( JH. GT. 176. AND. JH. LT. 186) 179,188 

—179 — rftGgJt  Jw-ire)  *10-  > tjt-  irs) 

IF(IA6E.GT.0.AND.IAGE.LT.100)500,188 



182  IAGE=  JH-176 

l£«ACt.-Dt-.41,ANO»IJU;f^I. 18-  1584,188  - ^ . 

188  IAGE=0 

PRINT  1000 

GO  TO  900 

.189  IF(IRECTYP.EC.31)  190,430 

190  lNhL!TlllT'*T?»»4)  VlNOLOdG) 

NUMVEC=(INOLD(18).AND.O70)*(2*»8)  + INO LO (19 ) » ( 2»* 4)  ♦ INOLO(20) 

1 XST  INOLO  < 2l  1-U87J481  -*  IJ4040132)  Y<  2«A81  * INOL«(8-31  • (2*  *41 

1 + INOLD  (24)) 

lYST  = (INOLO (25) *(2**12)  ♦ I NOL D ( 26 ) * ( 2**8)  ♦ INOLD (27) *( 2** 4 ) 

1 * INCL0(2e)) 

L=NOFTREC-IBLKNO 

430  ^?(  I^E^?YP,E  Q.^O  ) 435  ,4 5 0 

435  PRINT^436 


>*(?*»i.)*INOLO( 

$24) 

NUMBLKS=INOLO( 25) *( 2* * 12) ♦INOLD ( 26) *(2**8)*IN0LD(27)*(2**4)*INCLD( 

tZVi 

GO  TO  500 



495  IEOF=l 

PRINT-496 - 

500  RETURN 
END 


subroutine  UNFACRA  (JGGC.NWOtOIHfNSION  1080(89)  

COMMON  IBLKNO,ISTVEC,IAGE,INOLO(108Q) ,IRECRlS,L,NOFTR£C 
COMMON  XP,YP,RCP  (50 ) , SLOMER (2  0) .A(30),B(30)  .A&(10)  ,BB(10),  C(30), 

2 — DT70TTTT4B?T,Tr5D52“,T'5254,DEGLDNGr50)  ,DEGLAT{5D) 

INX=0$OO  2 J=1,NW0$00  2 1= 1 , 15$ INX= INX+ 1 

8 il«)UU  ll(Xl-l,78,A,SHIFTtJ(;8C{8)  , -4*<  15-1)  ) 

10  FORMAT  (5022) 

9 FORMAT  (*C  INOLD*) 


ENO 


6999 

MC6G^T16  0,CM145  00  0,TP2. 

FTTJ,  fC7A,T> , 4 9999 , T . 

MAP(ON) 

VSN,TA££4^S81J9- 
RE QUEST, TAPE1,HY,L . 
VSN,TAPE4=03  0172. 


2096 


MCCOMISH 
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0.  DKM  (II  l>\UI  l\ 


(l.l  0\<T\i<-» 

The  purpose  of  PAHT  IV  of  the  LORAN  Grid  Prediction  Program  is  to  update 
the  input  data  base  to  PART  II.  By  comparing  the  time  differences  at  measured 
and  calculated  fixed  points  in  the  service  area,  corrections  are  applied  to  the  pos- 
tulated top  soil  conductivity  values  to  zero  out  the  difference  between  the  former. 

A technique  described  by  Elkins^^  has  been  utilized  which  employs  a continuous 
Kalman  filter  with  updating  of  references  state  vector.  The  details  of  this  program 
are  described  in  Program  E.STIMAE.  It  is  pointed  out  that  of  the  nine  param- 
eters required  to  calculate  the  surface  impedance,  only  the  two  predominant  top 
soil  conductivities  are  altered. 

The  approach  taken  for  system  tuning  or  data  base  updating  is  as  follows: 

(1)  Compute  time  difference  in  PART  111  of  Program  for  locations  for  which 
measured  data  is  available. 

(2)  Identify  two  values  for  top  soil  conductivity  that  occur  most  frequently  in 
the  geographical  area  of  interest. 

(3)  The  variation  in  the  LORAN  time  difference  due  to  a variation  in  topsoil 
conductivity  is  computed  by  changing  in  turn,  each  of  the  two  conductivity  values 
where  they  occur  in  the  geographical  area  of  interest  and  proceeding  with  the  step 
outlined  in  " 1"  above. 

(4)  The  difference  between  the  observed  and  calculated  data,  the  rates  of 

change  in  time  difference  due  to  changes  in  each  conductivity  value,  estimate  of  the 
error  in  the  initially  chosen  values  of  conductivity,  and  estimates  of  the  error  in 
the  observations  are  input  to  a Kalman  filter  program  to  produce  an  improved  esti- 
mate of  the  value  of  conductivitv.  i 

I 

(5)  The  improved  values  of  conductivity  are  written  into  the  data  base.  ji 

For  the  test  program,  the  conductivity  was  changed  by  10  percent  to  form  the 

required  partial  derivatives;  absolute  conductivity  was  estimated  within  a factor  of 
two,  and  the  standard  deviation  of  the  time  difference  measurement  was  assumed 
to  be  100  psec.  Utilizing  forty  test  points,  this  technique  reduced  the  prediction 
errors  by  33  percent. 

A listing  of  the  program  is  presented  below. 


22.  Elkins,  T.  E.  (1976)  Empirical  Correction  of  Soil  Conductivity  Model, 

RADC/ET  Private  Communication. 

23.  Program  ESTIMAE  - Problem  No.  4723  (1975)  Analyses  and  .Simulation 

Branch,  AFCRL  dated  1 November  1975. 
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6.2  Pruxrum  Li^tin^  fur  Part  l\ 

I  ESTMATE.  7<»/7‘*  3PT  = L FTN  Ofc/io/77 

PRO&kA  6 ESTMAIE  <IH^JT=‘*31B,  OJIPU  1=1,0  13  I r APE1  = <, 318,  r APE  2 = 4013, 

1  r APE3=4tia,r APE4=*  oi3,TAPEio, r APEi i,r  APE2o=ii3ia) 

C....  UIl/EN  (1)  A bEI  OF  DIFFERENTIAL  EQUATIONS  OEbCRIBING  A OYNAHICAL 
C SrSIEn,  (2)  AN  ESTIHAfE  OF  I HE  STATE  OF  fHE  SYSTEH  AT  SOME  INITIAL 
C TIME  ALONG  HITH  STATISTICS  PROVIDING  A MEASURE  OF  THIS  cSTIMATt, 

C AND  (3)  A OATA  FILE  CONTAINING  OBSERVATIONS  BEARING  ON  THE  STATE  OF 
C THE  SYSTEM  AT  SPECIFIED  POINTS  OF  TIME.  COMPoTc  THE  BEST  ESTIMATE  OF 

C THl  state  of  the  SYSTEM  ALONG  HiTH  SOME  STATISTICS  GIVING  A MEASURt 

C OF  THE  ESTIMATE  AT  EACH  OF  THE  SPcCIFIEO  POINTS  OF  TIME.  THE  BtST 

C ESTIMATE  IS  BASED  ON  ILL  OBSERVATIONS  PROCESSED  UP  TO  THE  CURRENT 

C TIME.  A BACKWARD  SMOOTHING  -EATURE  IS  PROVIOEG  AS  AN  OPTION  SO  THAT 
C THE  ESTIMATlS  at  THE  ^IRST  < POINTS  OF  TIME  WILL  REFLECT  THE 
C INCLUSION  OF  ALL  OBSERVATIONS  UP  TO  THE  < TH  POINT  OF  TIME. 

common/case/case 

COMMON/RK/T  I,  T F,MR«l,S78,  SOT  ,DT,TOL,  SP,  NERRl,NERR2,ORO, 

I NSTEP.NKEJ ,ISJ 21 
COMMON/ST REF7 IXR,X  R(  1) 

COMMON/STESTM/IXE,  <E(1) 

COMMON/STRSOL/IX,X (1 ) 

COMrtON/STCOV7IP, JP,P(1) 

COMMON/STNCOV/  IQ.OCOV 
COMMON/OBSH/IH,  JH,  HI  1 ) 

COMMON/OBSY/IY,Y(1I 
C0nM0N/0BSC0V7Ik,R 1 1 ) 

COMMON/MEASCOV/RCOV (6) 

COnriOM/PARAMS/KSTRSDL  .HSTEP.TImETOl  ,K»  JNCHjKREH  IND.TOLRNCE,  lOER, 

1 KPLOT ,KENO ,<ONE ,<S , KSMUDT H, MA XS , MA XtC  , MPTY , N GPS ,NL RSR , NLRS W, 

2 NULF.NFILB.NSAT,  H N S AT , N GPSl , NGPS2  , NR  S E Q , NRR  AN , NWS  E Q , NWR  A N,  NPmAN, 

3 lOB, JZH, 1ST, I SI , I S2 ,KCRIT,KCBEGIN, KPR I N T ,KHR IT E , KO AT  A , KT Y? E S, 

. MAXO ATA,I A ,JA, KA, IM, JM, <M, MO, KPL WD,ir  Y PES(2)  ,PRINMOM( 3)  , KQKLOCK, 

5 lUERSET 

C....  THc  FOLLOWING  COMMON  AND  EOblVALENCE  STATEMENTS  ARE  DEPENDENT  ON 
C THE  DATA  FILE  USED. 

COnMON/BUFF/.BUFFI  .)  ,TIME 

DIMENSION  IDIl)  S EQUIVALlNCE  (D,  ID) 

COMMON  0(36),XRil(l>)  ,S(l->4) 

DIMENSION  Xa(12) ,DIAGJ (12) 
tXTERNAL  OEk 

NAMED  1ST  / VALUES /CA  SE  ,TSHRI  , T ST  OP  ,X  u , 0 I A G 0 , KPLOT,  KPRINT, 

1 <WR1Tl,KKE  wind, KSiUOTH, HAXS, KCRIT, MAX  KC , KCBEGIN,  KT Y PES , M AXD AT  A , 

2 lA, JA ,KA, IM, JM,<M,  NRK45/8,STEP,TOLRNC  E,SP,NtRRl,NERR2,0R0,  1ST , 

3 I El , IS2  , JZH,KSTRSOL, QCOV,RCOV,T IMETO.  , KP ARTS ,KPL HD , MO , RUNMA X, 

. KQKLUOK 

C...,  kEAU  in  INITIAL  PARAMETERS.  CALL  SUBROUTINES  TO  SET  INITIAL 
C PARAMtTERS  AND  TO  SELECT  UPIIONS. 

1 CASE=CASE+1.  I READ  VA.UES 

ificasc.lt.l.  ) Hb.a 

2 call  FIlEIOID 
CALL  observe 
CALL  STATS 

CALL  PACKIKPLWD, I,  1) 

SUT=ST£P  $ TOL=rO.RNCE 
IDtR=0 

:....  SEl  REFERENCE  jTATE  VECI3R  AND  DIAGONAL  STATl  COVARIANCE  MATRIX  TO 
C INITIAL  CONDITIONS. 

CALL  OIAGNALIIST ,X} , XR,0I  AGO,  IP,  P) 

C....  SET  MAXIMUM  RUN-TIME  LIMIT. 
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CL0CKQ=i;EC0N3(T) 

:....  Pi<iNT  initial  PAKAILTcKS, 

P«1NT  VALUES 

C>«a«  :}c.I  START  ANJ  sTUP  TINESi 
1 I=TSI ART 

i5Ji=isT*i  t isr ’1=ISI* ISfl 
C....  RiAU  iN  THt.  VEXT  DATA  KE:0RD. 

3 CALL  EILEIOKTI),  T E TURN  i { 4 , -.2 ) 

4 IF ai Mc.GT. TST0P)5, 8 
3 REND=1  $ G3  TO  -*2 

C....  TihE  OF  DATA  KECOP.J  LIES  WITHIN  RtQUESTEO  INTERVAL. 

6 I=IIM£  $ K3P=0 
il  TF^T 

HSiEP  = Aes(TF-Tn 

C....  IcST  MHt.TH£R  OR  NOT  TI  A40  IF  ARE  A PPR  0 X I MATiLT  EQUAL. 

14  IF (HSTtP.Lc.riHETO.) 21il> 

C....  iNTcGRATE  FGRHARO  TO  TF. 

15  IF<K0P)1'3,1c 

C....  map  THl  state  COVARIANCE  MATRIX  FORHARO. 

It  call  MATRIX (0, 1ST, IST ,0,’ , ip, XR( ISI  1),  1ST ,0,0 ) 
lbtkS£T=2  $ CAL.  RK4523 (XF ,DtR) 

CALL  MATRIX (0, 1ST, [ST,u,XR(lSTl)  , 1ST ,3 , I P ,L , j ) 

GJ  TO  21 

19  IF  (KRcSt-T.  tQ.  1)  GJ  TO  20  $ KR£SET  = 1 

CALL  MCVLt V (XR,XR.  , 1ST) 

20  call  RK4576 (XR,DE<I 

21  II=TF 

c....  select  the  uata  from  the  RECORC. 

IF (KOP.lQ.w )22,2S 

22  call  OHSEkVI  i I33=J  S I T Y PE S (1 ) =l  I T ° £S ( 2 ) = IH 

23  KOP= KOP*l  t iF (<3P.GT.<0ATA)36,24 

2<.  l.  = nU*  (KOP-1  )H  S T = 0(N)  S 1 = ID(N*2) 

IF(KOP.LQ.i)3u,ll 

C....  INITIALIZE  STATE  TRANSITION  MATRIX. 

80  UO  14  M=i,ISI  t .=IST*1  $ 00  9 K=l,IST  $ J-KpL 

9 XR(J)=4.  $ J=M+. 

iL  XR ( J ) =1. 

KR  lSl  T =0 

i.jLF,S£T  = l $ GO  T3  11 

C....  lNTER  APPROPRIATE  5U8KJOFINE  TO  FIND  THE  COMPUTED  OBSERVATION,  THE 
C OlFFiRENCt  eETwEEN  THE  GIVEN  OBSERVATION  AND  THE  COMPUTED 
C 03SERVATICN,  AND  THE  >ARTIA.3  OF  THE  COFPUTEO  03StRVAT 1 ONS. 

25  iF  (I.EQ.l)  26,27 

26  call  SStNSO^^KPARI S,D(N)  ,S)  $ GO  TO  35 


GO  TO  35 


26  call  SStNSO^^KPARI  S,D(N)  ,S)  $ GO  T 

27  iFd.  to. 2)20,29 

28  CALL  ESENSOR(<=>ARI  5,  D(N)  , S)  $ GO  TO  35 

29  1 F ( 4 . L E.  I.  A NO  a I . Ll.  ■ b ) 3 .1,  5 1 

30  call  MAGDA) A(KPARI>,D(N)  ,S)  $ GO  T 

31  iF (I. EQ.lb) 34.3b 

34  CALL  ACCELIKPARTS,  ) (N) ,SI 
....  SAVE  uESlRVATIjN  types  FOR  PRINTOUT. 

35  IF (KPRINT. and.  312)  5 30 ,23 

534  J=3*K0P-1  $ iNC33E (lb,>31,S)  J 

531  FORMAT  (iLH(Tl,tAl^,T,lL,.H, 12)) 

lNCOUE  (24 ,S  ,1TTPES)  ITYPES,! 

G3  TO  23 

30  1 I=TIME  f IF(IOH) 37,41 


73 


I ibriat i 


7^/ 


DPT  = L 


FTN  ‘t.S+i.li. 


37  IF  IKktSt-T.tQ.l)  C&.L  i107.  E J ( XRl, , Xk,  ISI  > 

3....  lNUK  the  FIlTER. 

IFl<0N£)38,39,39 

38  NFIL3=NFILa+l  $ iO  TO  tO 

39  NFILF=NFILF+1 

90  CALL  filter,  R£TUR'4S(99) 

C....  COMPUTE  THE  3EST  EiTIMATi. 

CALL  MATRIX(21,I3l,l,J,XR,IXR,X,rx,XE, IXE) 

C....  DlIcKMINE  PklHClPA.  /ALUIS  OF  ANGLcS. 

9000  bO  ^00  1=1,3 
90J  Xi  (I)  =PF.1NVA.  (Xl(1»  ) 

C....  REINITIALIZE  SiATE  REFERENCE  VECTOR. 

CALL  MCVLEV (XE,XR, 1ST) 
c....  chlCk  for  maximum  ?un  time. 

CLbCK  = St.COND(T>-C.3CK0  5 I F ( CL  OCK  . Gf  . RU  NMA  X ) KENO  = l 
C....  GAINER  APPROPRIATE  STATISTICS  ON  THc  ESTIMATION. 

call  STATSl 

IF (KENC. Ed. 0)  CALL  FILtI32(TI),  RET  URN S ( 9 , 92) 

C....  initiate  end  PROCESSING.' 

92  CALL  FILEI03(T) 

93  GO  TO  1 

99  PRINT  95,NGPS,KUNE 

95  FORMAT  (*0SINSULAR  MATRIX  E NC  OUNT  £ RE  (J  .»  5 X • GRO  J P = * I 9,  6X  • KONE=»  12 ) 

SIMPLY  cQUATE  BtST  ESTIMlTt  VECTOR  ANO  STATE  REFtkLF.CE  VECTOR. 

91  call  MOVLlV (X<, XE, [ 51 ) i 00  950  <=1,I3I 

950  X<K)=l.  $ SO  TO  ,000 

9d  CONTINUl  $ END 


N PkINVAl  79/79  gPT=1  FTN  9.59.,m 


FJNCTION  PR  INVAL  (ANGLE) 

C....  PkINVAL  K£TU\NS  THl  PkInSIPAL  VALUl  0-  The  given  angle. 
C ( - “I  .LT.  PnINVAL  ..E.  pi  ) 

DATA  Pl/3.A9l5a2b3358979J/ 

A= angle 

1 AAiAbS(A)  $ iF(A».LL.Pt)  GO  TO  2 
A= A-SIGN(2. *PI ,A)  $ GO  TO  1 

L PRINVAL=A  S ivcTJRN  S END 
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SUaKOUTlNE  SrAIS 

....  RECORDS  STATISlICi  IN  ORJER  10  Ei/ALUAIE  THE  PERFORMANCE  OF  THE 
ESIIMATIOK  PROCESS.  :ONTKO.S  THl  INlIIAriON  OF  THE  BACKWARD 
SMOOtHiNG  PROClSS  ANJ  THE  ASSOCIATED  RANiOM  FILE  WRITING.  CONTROLS 
ALL  PUNTED  OUTPUT  ON  FILcS  1,E,3,  ANO  4. 

COMMON/STREF/IXR,XR(1) 

COMMON/SlESIM/IXc.,  <L  (1) 

COMHON/STRSDL/IX,X ( 1) 

COMMON/STCOV7IP, JP, P (1) 

COMMON/SICOWO/IPO, D (1) 

C0MMON/08SY/IY,Y  (II 

COHMON/PARAMS/IQ (7  I ,KPLOI  , KE NO , K ONE , K> , <SMOOTH,MAXS .MAXKCiMPTY  , 

1 NGPS ,NLkSR,NL<SW,  NF ILF , N FI L B , NS  AT, NN> A I , NGPSl, NGPS 2 , NRSE 0, NRR AN , 

2 NWSEQ.NWRAN, NPRAN,I0B,JZH,ISI ,IS1, IS’  , < C R IT , KC BE G I N , KPR INT , KW R I T E 

3 , Ji3(  11)  .OBTYPES  (2  I 
COMMON/CRIT/CRIT (II 

C....  T Ht  FOLLOWING  COMMON  ANO  EQUIVALENCE  STATEMENTS  ARE  DEPENDENT  CN 
C THE  DATA  FILE  USEO. 

CDMH0N/6UFF7LdUFF(.) ,M(1’7) 

dimension  XT(1)  $ EQUIVALENCE  ( W, T IW E ) , (W ( 1 Q> , X T I 

CO MM  UN  DX(12> ,X0 (12) ,IS(24) ,0  (50  ) 

DIMENSION  S(l)  i EQUIVALENCE  (IS,S) 

LOGICAL  NRSEQ,NRRAV, NWSEQ, NWRAN, NPRAN 
DAIA  MAXK75 7, MORE/ iH  M0vE7 
C....  SET  INITIAL  CONSTANTS. 

NLRSR=NLRSW=NGPS=N- ILF=N- ILB=NSAT=NNSAI =KST=KEND=C 
NRSi.Q  = .T.  $ NkRAN=NWSEa  = NWRAN=NPRAN=.F. 

KS=;KSMOOTH  $ KONEsl 

C....  SET  PRINT  OPTIONS.  AT  MjST  1 OPTION  ’ER  FILE  IS  ALLOWED. 

C KPkiNT  BIT  0.  ORiNT  TIME,  STATE  VECTOR  ON  FILE  1. 

C KPRINT  BIT  1.  PR. NT  fIME,  STATE  kcSlUUA.S  VECTOR  ON  FILE  2. 

C KPRINT  BIT  2.  PRINT  IIME,  STATE  DEVIATIONS  VECTOR  ON  FILE  2. 

C KPRINT  BIT  3.  PRINT  TIME,  STATE  DERIVATIVES  VECTOR  Ou  FILE  2. 

C KPRINT  BIT  4.  “RINT  TIME,  JIAGONALS  OF  STATE  COVARIANCE  MATRIX  ON 
C fill  3. 

C KPRINT  BIT  5.  =RINT  TIME,  STATE  COVARIANCE  MATRIX  ON  FILE  3. 

C KPRINT  BIT  b.  PRINT  TIMl,  fORQUES  ON  FI.E  3. 

C KPRxNT  BIT  7.  PRINT  COMPLEIE  STATISTICS  ON  FILE  •♦. 

C KPRINT  BIT  e.  PRINT  ’ARTIAL  STATISTICS  ON  FILE  4.  (KCRIT.NE.L) 

C KPRINT  BIT  9.  ’RINT  TIME,  SOME  STATISTICS,  ANO  OBSERVATION  RESIDUALS 

C VECTOR  ON  FILE  4. 

IF (KPRINT) 2 ,1 

1 ASSIGN  56  TO  Ml  $ GO  TO  17 

2 CALL  WKITER0(4,0,3I 

ASSIGN  41  TO  Ml  $ ASSIGN  46  TO  M2 
ASSIGN  52  TO  M3  I ASSIGN  55  TO  M4 
ASSIGN  641  TO  M41 

C....  SET  PRINT  LIMITS  ON  STATE  VECTOR  ELEMENTS. 

I52=MINl (12,IS2, IsT , (IX-ISlNl) ) $ IS  4= I Sl  + I S ’-1  I IS3  = IS2tl 

PRINT  16, (I  ,I  = IS1, I S4) 

PRINT  160  $ PRINT  161 

1 = J=1 

3 K=KPRINT.AND. J $ IFCK.NE.O)  GO  TO  ( 5 , 5 , 7 , 9 , Id , 11 , 1 3 , 14 , 14, 15 ) , I 

4 Jj2*J  $ i=l41  S IFd.GT.lO)  17,3 

5 ASSIGN  40  TO  Ml  S GO  Tj  4 

6 ASSIGN  42  TO  M2  S GO  TO  4 

7 ASSIGN  43  TO  M2  t GO  TO  4 
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■5  ASblGN  1,4  TU  M2  { GO  TJ  4 

10  AbSlGN  47  TO  M3  S GO  TO  t 

11  AiGIGN  Hb  TO  M3  S CO  Tj  4 

13  AbSIGN  51  TO  M3  5 GO  TO  4 

14  KSi=l 

l40  ASbIGN  54  TO  Mf 

iF(i.LU.G)  Assign  33  to  m-*  s go  to  « 

15  ASSIGN  544  TO  M41  $ GO  TO  l4i, 

lb  FORMAT (*lO£SCRiPTI On  OF  >RlNTOOTS  ON  -l.ES  1,  2,  OR  3»/5X*ONLf  IriE 

1 FOLLOHING  COMPONENTS  OF  THE  STA1E  VEOTOR,  STATE  RESIDUALS  i/ECTOR, 

2 sIATt.  DE7iATIuNS  i/ECTOk,  STATE  UEr1i/ATI7ES  7ECTOk  Ok  TH£*7»  FOLLO 
3H1NG  COLUMNS  OF  THE  STATi  COi/AkIANCE  IATkIX  OR  THE  FOLLOWING  DIAGO 
ANAL  ELcMcNTS  OF  THE  STATi  COVARIANCE  NAIRIX  ARE  PRINTED  0U1.*/1X, 
51215) 

ibi,  FORMAT!*  FOR  IHl  jiSERVAilON  RESIDUALS  VECTOR,  ALL  CCMPONENTS,  AS 
IWiLL  AS  THt.  OBslRVUION  tYPtS,  ARE  PRINIlO  OUT.*) 

161  FORMAT  I 

a'-UEsCRIPTION  OF  jIATISIICS  PRINTOUT  ON  FILE  4* /5X' GP . . . . POS IT  ION 
7UF  uATA  GROUP  WRT  -IRST  JATA  GROUP  PRO C E S SED* /5 X* LR  R . . . . L uC AT  ION 
aOF  logical  RECOkD  OF  UAIi  WRT  FIRST  LOGICAL  RtCORD  READ*/5X 
3*LR  W....CUMULATIV1  NUnBER  of  LOGICAL  RECORDS  WRITTEN  TO  OUTPUT  DA 
ATA  FILE  */3X»F  F . . . . CUMU L AT  I VE  NUMBER  OF  TIMES  ENTE»EO  FIL 

GlcR  FOk  FORWARu  Fi.TEI<iN3*/5X»F  B.,  . .0  UMULATI  Vl  NUM6ER  OF  TIMES  EN 
CTcklU  filter  for  dlCKWAR^  SMOOT H IN3 *73 X * S AT . . . . CU MU L A T I VE  NUMBER  0 
OF  TIMlS  PcRFORnANCE  CRIT.RIA  S AT ISF I EO * 7 a X*NS AT  . . . . C UMUL A TI V E RUMB 
EER  OF  TIMES  PERFOR1ANCE  ^RITERIA  NOT  > A T I S FIE  0* 7 5 X » E LE ME  NT . . . , R t L A 
FIIVl  location  within  STAIE  OR  STATl  RESIDUAL  V_C  I OK* 7 5 X* V AL U E. . . . T 
GlSI  VALUt  which  EXOlEDEO  SPlCIFIEO  CRirERIA*75X*MORE....FLAG  IMPLY 
HiNG  THAT  MORE  -LEHENTS  FAlLtO  TE ST* 75< » T I ME. . . . Tl ME  OF  DATA  GROUP* 
I) 

17  Assign  ie  to  Mo 
C....  SET  CRITtRIA  OPTIONS. 

C KCrIT  bit  c,  covariance  lihit  criteria. 

C KCRII  BIT  1.  COMPARISON  WI 1 ri  EXISTING  STATE  CRITERIA. 

ASSIGN  37  TO  M5 

iF (KCKIT.lQ.I)  assign  21  TO  M5 
IF IKCRIT.EQ.2)  assign  22  TO  M5 
C....StI  SEQUENTIAL  WRITE  OPTIuN. 

IF (KWRIT t.EQ. 1 ) NWSLQ=.NjT.NWSLQ 
rET  URN 
ENTRY  SIATSI 

0....  Usual  entry  POINT  INTO  STATS. 

00  TO  Mb, (16, 2u) 

C....  WHlN  PkOCESSINj  F.RST  OAIA  group,  specify  necessary  OPTIONS. 

16  assign  lL  to  Mo 

IF  (KPLOT.EQ.1)  NPRAN=.N jT .NPRAN 
IF (KSMOOTH) 19, LU 

Ij  NWkAN=.N0T. NWRAN  E NWS . 0 = . N OT . N WS EQ  E NPRAN  = . NOT.  NPRAN 
LGPNS=G 

C....  CHECK  WHETHER  DR  NDT  CRITERIA  ARE  SATISFIED, 

2l  K=l  t DO  2b  1=1, 1ST  E CO  TO  M5, (21 , 22 , 37) 

21  J=1*1P*(I-1)  $ VAU=P(J)  f GO  TO  420 

22  vAL=Xt  (I)-XT( I) 

220  IF (A8S  (VAL)  .Lt. CRir ( I) )2o, 23 

23  K=K*1  S IF(KSI)2.,36 

C...  sAVl  data  on  LLlMENTS  not  SATISFYING  CRITERIA. 

24  IS(K)=I  S K=K41  E S(<)=VAL  $ IF( K , G T. 2* MAXK) 2 5 , 2 E 
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2S  K = K-1  1 IS(tC)=MLI?E  S GO  ro  36 

2c  COi.IiNUL  $ IF(<)  56,27 

27  NbAIiN6AT.l  $ I? ( KS. £Q. 1 ) c8, 39 
2 8 IF  ( INGPS-L  GP'4S)  .Gc.  rtAXKCl  29,31 
G....  lAXKC  CONScCJTitfc  > A TI SF  A C T Okr  ESUMAIlS. 

29  IF iLGPNSlou  ,35 

C....  oAcK.<ARO  bMOOIilNS  NECESSARY. 

3u  Ni,P31  = NGPS-KG8cGIY  $ N.  RSR=  NL RSR- KC3  E GI  N $ GO  10  33 

31  IF .NGPS.Gc. MAXS) 3c, 39 

32  N3Pil=NGPS 

33  X0Nc=*1 

3.  NGPS2=NGPS  t KS--1 

..R8tJ=.N0T.NRSEQ  5 iYRRA  N=.  NOT.  NK.RAK  S GO  TO  39 
B-OKHARL)  8MOOTH1U5  NOT  NECESSARY. 

35  N&PS1  = 1 $ NLRSR='ILRSR-'IGPS.1  I GO  10  3A 

C....  AN  ONSATISFACToRY  ESTIMATE. 

30  NNiAl =NNSAI +1  $ .&PNS=NGPS 

37  1 F (KS. EO. 1) 3l , 39 

lXIRACT  OIAGONAL  POrxIION  OF  COVARIANCE  MATRIX. 

39  00  39b  1=1, IS T 
39j  0 ( I) =P  (I  + IP*(I-1  )) 

gale  MATKlX(22,lSr,l,iJ,XE,iST,XT,IST,<D,IST> 
c....  PRINT  Selected  in-ormation  on  output  -iles. 
b(l)=TlME  S GO  ro  Ml, (.0,-1, 55) 

AL  GALL  MCVLtV (XE (ISit  ,0(2) , IS2)  $ CAL.  HRIT ER (1 , 0 , 1 S3 ) 

A1  GO  TO  Me  , ( A E,  4 3,  AA , A 6 ) 

a2  gall  MCVLEV  (X(ISl)  ,0(2),IS2)  $ 60  TO  ts 

A3  CALL  MCVLLV(X0(IS1) ,0(2) , ISe)  $ GO  TO  a5 
AA  GALL  MCVLEV  (OXdSl)  ,0(2)  ,IS2) 
a5  call  mRITER (2, 0, 135) 

AO  GO  TO  M3,  ( A 7,  Ac  , 51 , 520 

A7  CALL  MCVLcV(G(1S1) ,0(2),ISe)  $ L=IS5  $ GO  TO  53 

A6  L=l  $ 00  A9  I=I51,ISa  $ DO  A9  J=1S1,ISA  t L=la1 

A9  0(L)=P(iAlP*(J-i))  t GO  TO  5t 
5:  CALL  WRITER  (3,0, L) 

51  CONIINUE 

52  GO  TO  Ma  , ( 5 3,  5a , 55  ) 

53  . F (K ) 5a, 55 

5.  CALL  MCVLEV(NGPS,0(2) ,7)  $ GO  TO  MaI , ( 5 AO , 5 A1 ) 

5AL  K = l $ IFdoe.cQ.  0)  GO  10  5a2 

call  MCVLEV(08TYPE>,0(9)  , 2)  S CALL  1 0 VL E V C Y ,0 ( 1 1 ) , I 0 6) 
K=eAIOB  $ GO  TO  5A2 
5A1  call  M0VLEV(1S,U(9) ,K) 

5a2  call  WfinER(A, 0,8*0 

55  IF (NPRAN.ANO, .not. YWRAN) j6,57 

C....  SAVE  DATA  FOx.  tVENTJAL  P.OTIED  OUTPUT. 

56  GALL  MOVLeV  (XE,0(2( , 1ST)  $ I = 2aIST 
lALL  MCVLeV(X,0(I) ,IST)  $ I=I*IST 
CALL  MCVLEV (XC,0 (I ), aST)  $ I=I+IST 
CALL  MOVLEV(0,0(I),ISr)  $ I=I*IST-1 
CALL  PACKK  I.O.MPT  f ) 

57  RETURN  $ ENO 
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SUBROUTINE  FILElO(r),  REF  URNS (NR1«NR2I 
C....  THIS  ROUTINE  IS  TTE  CONTROLLING  ROUTINE  FOR  ALL  FILE  INPUT  AND 
C OUTPUT  IN  THE  PROGRAH. 

COMMON/STREF/IXR,XR(l) 

COHHON/STtSTH/IXEi  <E (1) 

COMHON/STCOV/IP, UP, P (1) 

COHMON/STCO VO/lPO,  3(1) 

COMMON/RANFILE/HAXiS,MS(l ) 

COHHON/PACKEO/NLRSi HHOROS,MH(5ie) 

E3UI VALENCE  ( HH, IH 4 1 ) i (MR ( 2) , IHH2) 

COMMON/PARAMS/KSTRJOL.HSr EP.TIMETOL ,K» UNC H, KRERIND,  L RMOCNT ( 2 ) , 

1 KPLOT,KENO,<ONE>K>  >KSH03TH,HAXS»HAXK: , NPTY,NGPS|NLRSR|NLRSW, 

2 NFILF,NFILB,NSAT,  4NSAT, NGPSL , NGPS2  >NR SE QiNRRAN, NMSEQ>  NMRAN, NPRAN, 

3 IOB,JZH,IST,Zl(22)  iKQKDOK 

C....  THE  FOLLOWING  CONHON  AND  EQUIVALENCE  STATEMENTS  ARE  DEPENDENT  ON 
C THE  DATA  FILE  USED. 

COMMON /BUFF /length I, LENGf  H W, M AXL RH, LRW , W(127) 

E(lUIVALENCE  (H.TIME) 

COMMON  RECORD(14) 

LOGICAL  NRSEQ,MRRAN,NHSE3.NWRAN, NPRAN 
DATA  KOPEN/-1/ 

C....  SET  FILE  INPUT/OUT>UT  OPTIONS. 

IF (KSMOOTH.NE. O.OR, XPLOT. NE.O)  KOPEN=<OP£N+1 
IF(KOPEN)101, 102,1)1 

102  CALL  OPENMS(20,MS, 1AXMS»1 , 0)  S KOPEV =<OP£N+ 1 
101  MAXS=MIN0 (MAXSiMAXiS)  S NAXKC= NINO (M AX KC,NA XS ) 

IF (KREHINO.NE. 0)  REMIND  10 
ASSIGN  36  TO  Ml 
KREWRIT=0 
RETURN 

ENTRY  FlLEIOl 
C....  RANDOM  FILE  READ. 
lOO  IF(NRRAN)1,15 

1 IF(KS>2,6,15 

2 IF(NGPSl.EQ.l) 3,6 

C....  BACKWARD  SMOOTHING  EITHER  NOT  NECESSARY  OR  NOW  COMPLETE. 

3 KONE=l  S KS=u  S ASSIGN  37  TO  HI 

NWSE 0=.NOT.NWSEQ  t NWRA N=.NOT. NWR AN  S NPRAN=. NOT. NPRAN 
IF(NGPS1.EQ.NGPS)37.1Q 
C....  CONTINUE  WITH  BACKWARD  SMOOTHING. 

6 IF(NGPS.EQ.NGPS1)3, 7 

7 ASSIGN  11  TO  Ml  $ GO  TO  10 

8 IF (NGPS.NE.NGPS2)  GO  TO  3 
C....  RESUME  SEQUENTIAL  -ILE  READ. 

NRRAN=.NOT.NRRAN  S NRSlQ  = . NOT. NRSEQ 
KRtHRIT=-l  B GO  TO  11 
C....  READ  IN  A LOGICAL  RECORD. 

9 NGPSl  = NGPSl+<ONt  t NLRjR  = NLRSRi-KONE 
13  CALL  REAOMS (20.LRM, MAXLRW, NGPSl) 

LRW=LRH-2*ISt 

NGPS=NGPS1  E GO  TO  Ml, ( 1 1 , 36, 3 7) 

C....  MOVE  APPROPRIATE  VALUES  FROM  W ARRAY  TO  STATE  AND  COVARIANCE 
C ARRAYS 

11  I:LRw«l  $ CALL  N3VLEV(W (I),XE,IST) 

I = I«^IST  B CALL  M0VLEV(4(I),0,IST> 

CALL  OIAGNAL( IST,X£ , XR,0, IP,P)  B T^TINE 
ASSIGN  3b  TO  Ml  B GO  TO  100 
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C....  SEQUENTIAL  FILE  RE^O. 

15  IF(NRSEQ)16f22 

C....  kEAU  in  a losical  record. 

16  KEAO(IC)  LRH,  (Ntl)  , I = 1,LRH)  $ NLK.SR=  N. KSk^l 

c....  crlck  last  binary  reao. 

IFdO  EOF  PE<1U)  J2  + , 17,13 

17  NLRSR=NLRSR-1 

22  KENU=1  $ RETURN  YR2 
C....  DATA  LOGICAL  RECOR). 

Zk  IF  (TIME.LT.T-riMEr  35 

C....  ACCEPTABLE  DATA  RiCORO, 

35  N3PS  = NGPS«-1 

36  RETURN  NRl 

iNTRY  FILEI02 

C....  RAUOOH  FILE  NRiTE. 

37  IF (NWRANlAfl  ,41 

40  1=lRH+1  $ CALL  NOVLEY(<E,W(I),ISr) 

I = itIST  $ CALL  r<37L£70,W(I)  ,IST) 

LRW=LK4+2*ISI 

CALL  WRITHS <20,LRW, LRH»1, NGPS,KREHRITI  5 GO  TO  49 
C....  SEQUENTIAL  FILE  HRCTE. 

41  IF (NWS£Q)42,45 

42  I=LENGTHR+1  $ CA.L  MOV. E V ( X E , H ( I) , IS T ) 

I=I*LtNGTHW  $ CA.L  MOV. E V (0, W < I ) , IS!  ) 

LRH=LENGTHR+2*LEN3r HW 

WRITE(ll)  LRH,  <W(It  , I = 1,.RW)  $ NLRS>(  = NLRSM«-1 

C....  OUTPUT  DATA  FOR  EVENTUAL  PLOTTING,  COMPARISON,  ETC. 

45  IF (NPRAN)4b,49 
4b  IF  (MPTY.EQ.l)  48,49 

48  M=NLRS*1 

CALL  NKITHS  (Ell  ,NM,  VWOROS,  M,KREWRIT)  $ NHbRDS=0 
IF (H. EQ.MAXMS)  N PR A N= . NOI  . NPRA N 

49  IF(KENO.EQ.1)22,1jJ 
cNIRY  F1LEI03 

c....  perform  end  processing. 

5u  F3RMAT (♦-MINIMUM  AVO  MAXIMUM  VALUES  F3R  OATA  ON  PLOT  F IL E . • T 7 0 • WO R 
105  PER  GROUP*I<„TliO*TOrAL  LOGICAL  RE:0R0S*I4/) 

51  FORMAT (I4,2X,1P2£23. 12) 

53  IF (KPL CT) 54,56 

54  IF(NHORDS.EQ.O)  GO  TO  55  $ M=NLRS+1 

CALL  HFITMS (2u , WW, V WOKUS , M , KRE MRI T) 

55  call  PACK2(N,W,MPrO 

PRINT  5C,IWWl,IHH2  S 03  55li  I=l,IWMl  $ M=I*2  5 N=M*IWW1 

55u  PRINT  51,1, WW(M) ,WV(N) 

call  H PITHS (20 ,HH, VWOROS.l ,KREWRIT) 

56  IF (KQKLOOK) 57,62 

57  UO  q1  1=1,4 

lNOFILE  I $ BACKSPACE  £ $ BACKSPACE  I 

KEAO(I,56)  RECORD 

58  FORMAT (1 3A10, A7) 

IFdO  EOF  PEd))  59, 61, 61 

59  Print  6u,i 

64  FORMAT (•-FILE*I2*  -AST  RECORD*/) 

PRINT  56, RECORD 

61  CONTINUE 

62  RETURN  $ END 
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N lOEOFPE  74/7I*  OPr=l 


FTM 


FUNCriON  10  EOF  PE(J» 

C....  rms  FUNCTION  CHE3<S  for  NORMALi  END-DF-FILE,  AND  PARITY-ERROR 
C CONDITIONS  ON  UNIT  J F3LL0M1<1G  A BINARY  READ. 

DATA  KOUNT/0/ 

K0UNT=K0UNT+1  _ 

C....  NORMAL. 

K=-l 

IF (EOF(J))10tll 
C....  END  OF  FILE. 

IG  K = u S PRINT  2i]  t GO  r 0 13 

11  IF(I0CHEC(J>)12,1>* 

C....  PARITY  ERROR  I N OR  AFTER  RECORD  INDICATED. 

12  K=1  $ PRINT  21 

13  PRINT  22,J,K0UNT 
1<*  10  EOF  PE  = K 

2S  FORMAT  (2<*MJ  EYD  OF  FILE  •*•••) 

21  FORMAT t25H0***»*  PARITY  ERROR  •••*») 

22  FORMAT!*  UNIT  » 1 2, I S X* • RECORD  *14) 

RETURN  $ END 


E OBSERVE  74/7<.  OPT^ll 


FTN  4.5*414 


SUBROUTINE  OBSERVE 

C....  IMIS  ROUTINE  SELECTS  THE  TYPES  OF  DATA  REQUESTED  FROM  THE  DATA 
C FILE  AND  PLACES  THEM  IN  AN  INCREASING  SEIUENCE  IN  STORAGE  FOR  USE  BY 
C PROGRAM  ESTHATE. 

COMMON/SENSORS/SEMAdS)  .:OILS(3,2) 

CONNON/OBSH/IH.  JH.  M ( 1 ) 

COMMON/PARAMS7IQ(3i) (KOAT A,KTYPES,MAX3 ATA.IA, JA,KA, IH, JH,KH,HD 
C....  THE  FOLLOMINS  COMMON  STATEMENT  IS  DEPENDENT  ON  THE  DATA  FILE  USED. 
C OMM ON/BUFF /L BUFFI f) ,H(1;7) 
dimension  data (1) , lOATA(l) 

EQUIVALENCE  (M (1 9) . NO)  , ( V ( 2C ) • OAT A, IDA T A I 
COHHON/TORQUES/COI . (3) >T9RQHAG(3) .OENCOEF 
COMMON  D(36>,T(20)  ,L0C(2J) ,MN(4) ,LN1(4) 

DIMENSION  LM2(>»)  ,LM3(4) 

DIMENSION  C0ILSP(3, 2) 

DATA  C0IL,C0ILSP/9'a./ 

DATA  KMASK/ 777 7 7 77 1777771 7777 7 48/ 

LM2(3)=LM2(2)=LM2<l|)=JN  S L M3  ( 3)  = L H3  ( 2 ) =LM3  ( 1 ) : KM 

LM2(4)=JA  $ LM3(*)=KA 

RETURN 

ENTRY  08SERV1 

C....  SELECT  DATA  FROM  DATA  FI.E. 

KDATA=K=MM(4) =MM(3I =MM{2) =MM(1)=0 
KCOUNT  = (i 

LM1C3)=LM1{2)=LM1(1)=IM  $ LM1(4)=IA 

20  IF(K.EQ.ND)  SO  TO  J1 

J=MD*K*1  $ K=K*1 

I = KTYP£=10ATA(  J*2)  $ I-d.GT.e)  1 = I.  A NO.KMASK 

IFd.AND.KTYPES)  21,20 

21  IF (KTYPE.ANO. 8) 22, 230 
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7<./74  3PT=ll 


FTN 


C..  ..  STORE  COIL  OATA. 

22  L=KrrPE-7  $ UO  2J  1=1,1 

23  COlLSP(I,L)=OATA(JH)*COlLS(I,L) 

CALL  MATRIX (21, 3 ,1,3 ,CUILSP,3,COILSP(* ) ,3 ,COIL,3)  $ GO  T3  2b 
23b  IF(KTYPE.ANO.Id) 23L,24 
C....  SAVE  THE  ONE-HALF-THO-V-SQOARED. 

231  OENCOEF  = DATA(  J«^3)  S GO  TO  20 
C.*..  SAVE  OTHER  TYPES  3-  DATA  FOR  PROCESSING. 

24  1F(KIYPE.LE.2)  GO  TO  26  S L =MOO (KTY^ E , 9 > -3 
MM(L)  = MM(L) +1  S IF(MM(.).£Q,L(11(L))25,20 

25  LM1(L)=LM1(L)«LM3(.)  $ IF (L Ml (L  ). GT.  L N2  (L 1 1 LN1(L)  = 0 

26  KDATA  = K0ATA4'1  S T ( KOATA)  =0ATA(  J)  S L0C(K0ATA)=K 

XCOUMT=KCOONr»l  S IF<Kr YPE.EQ. 1)  KC3 U NT =KCOUNT* 2 $ GO  TO  20 

C....  LIMIT  DATA  I-  NECESSARY, 

31  L=KCOUNT-MAXOATA  i IF(..GT.O)  KDAT A= K3 ATA-L 

C....  RETURN  WHENEVER  NJ  DATA  IS  SELECTED  OR  NO  ORDERING  OF  OATA  IS 
C NEEDED. 

IF (KDATA-1) 4U, 36,310 

C....  ORDER  THE  DATA  IN  IN  INCREASING  SEQUENCE  BY  TIME. 

313  03  37  I=2,KDATA  t J=I-L 

T£SI=T(I)  S IFTTEST.GE.TTJ))  CO  TO  37  $ lO=LOC(I)  $ L=I 

32  J=J-1  $ IF(J)33,14 

33  IF (TEST.LT.TC J)) 32, 34 

34  J=J+1 

35  LN1=L-1  $ T(L)=r(LHl)  t LOC ( L )=LO; ( LM 1) 

L = L-1  $ IF(L.  EQ.  J)  36,35 

36  T ( J) =TEST  $ .0C( J»  =L0 

37  CONTINUE 

C....  STORE  THE  SELECTED  OATA, 

38  DO  39  1=1, KOATA  i L=MD* ( LOC ( I) -1) ♦!  $ J=H0*<I-1)»1 

39  CALL  HOVLEV(DATA(t)  ,0(J),l«0) 

40  RETURN  S END 


81 
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aPT::l 


FTN  *..5+41i* 


iUbKOUTINE  PACK(N,/,M) 

C....  THiS  kOUTINt  SrOk£S  DATA  BEFOftE  IT  IS  OUTPUT  AS  A SINGLE  LOGICAL 
C RECORD.  THE  RANGE  OF  THE  DATA  IS  ALSO  C3H»UTEO. 

DINE  NS  I ON  ;(!),/ HI '4(5iJ)>i/NAX(5G) 

COMHON/PACKEO/NLRS.NHOROS, VV(5 12) 

Eaui VALENCE  (VV(1) ,IWO),( VV(2)  ,J) 

C....  MAXIMUM  NUMBER  OF  VOROS  IN  A LOGICAL  RECORD  IS  512. 

c....  maximum  number  of  roros  in  a group  is  53. 

DATA  MAXO/512/,MAXV0/50/ 

C....  INITIALIZE  PARAMETERS. 

NLRS=D 

iH[)=N=MlNL  (N.HAXRJI  J ^ GP=  ( MAXO-2  ) /I  HO 
DO  1 K=1,IHD  $ V1IN(K) :l.t99 

1 VMAX ( K)=-l. E39 
NMOROS=J=0 
KEIUKN 

ENTRY  PACKl 
C....  STORE  THc.  DATA. 

IFIM.eQ.I)  j=o 

IF(J.NE.O)  GO  TO  3 $ M:  0 $ NLRS  = N.RS  + 1 

3 U=U«-1  $ N=MIN0  (N,  IWO)  S L = IHD*<U-t) +2 

00  2 K=1,N 

Vl=V(K)  $ IF(VT. LT.VMIy(K))  VMIN(K)=VT 
IF (VT.GT.VMAX(K))  VHAX(KP=VT 

2 VV(L<-K)=VT 
IF(U.EQ.NGP)  M=l 
NW0RDS=L+IW0 
RETURN 

ENTRY  PACK2 

C....  READ  OUT  THE  SIORiO  RANGE  VALUES. 

NHGR0S=2*IHD+2  $ U=NLRS 

CALL  MCVLEV (VHIN.VV (3) ,IhO)  S CALL  1 0 VL E V ( VNA X. VV ( 3 + IHO) « I WO ) 
RETURN  $ ENO 


E DIAGNAL  rH/ZA  OPT=L  FTN  A.5+A1A 


subroutine  oiagna. (n,X|X<,o,ip,p) 

C....  utility  subroutine  to  move  THE  N-DIMENSIONAL  STATE  VECTOR  X INTO 
C XX  AND  TO  MAKE  THE  N- 0 1 ME NS I ONAL  VECTOR  0 THE  DIAGONAL  MATRIX  P. 
C....  NOTo.  IP  IS  COLUMN  SIZE  OF  MATRIX  P. 

DIMENSION  X (1) , XX( I)  ,0  (1>  ,P(1) 

CALL  MCVLEVIX.XX.NI 


00  2 1 = 1, N 

$ L=I’» 

00  1 K=1,N 

$ U=<>L 

1 P(U)  = (;. 

U = ItL 

2 P(U)  = D (I) 

RETURN  $ 

END 
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